home *** CD-ROM | disk | FTP | other *** search
/ SGI Developer Toolbox 6.1 / SGI Developer Toolbox 6.1 - Disc 4.iso / public / SciAn / src / ScianSIMSFiles.c < prev    next >
C/C++ Source or Header  |  1994-08-01  |  100KB  |  3,397 lines
  1.  /*ScianSIMSFiles.c
  2.     File readers for Steacie Institute of Molecular Sciences
  3.     National Research Council of Canada
  4.  
  5.     Christina Lam 
  6.     Summer Student
  7.     August 1993
  8.     email: lamchri@ecf.utoronto.ca (until May 1, 1994)
  9.  
  10.     XYZ, GRD, GSN1, GSN2, GSN file readers
  11.  
  12.     XYZ: ball/stick molecules in XMol's XYZ file format
  13.     GRD: reads BioSym's DMol PLOT 3d grid files
  14.     GSN1,GSN2: molecular orbitals from gaussian primitives
  15.     GSN: takes Gaussian92 output, generates molecular orbitals
  16.          for standard N-M1G basis sets, generates ball/stick
  17.          molecule
  18. */
  19.  
  20.  
  21. #include "Scian.h"
  22. #include "ScianTypes.h"
  23. #include "ScianWindows.h"
  24. #include "ScianControls.h"
  25. #include "ScianTextBoxes.h"
  26. #include "ScianFiles.h"
  27. #include "ScianLists.h"   
  28. #include "ScianPictures.h"
  29. #include "ScianErrors.h"
  30. #include "ScianTimers.h"
  31. #include "ScianDatasets.h"
  32. #include "ScianIDs.h"
  33. #include "ScianStyle.h"
  34. #include "ScianSciences.h"
  35. #include "ScianArrays.h"
  36. #include "ScianTemplates.h"
  37. #include "ScianTitleBoxes.h"
  38. #include "ScianButtons.h"
  39. #include "ScianColors.h"
  40. #include "ScianSIMSFiles.h"
  41.  
  42.  
  43.  
  44. extern ObjPtr visBalls, visGeometryClass;
  45.  
  46. /* ---------------------------------
  47.     XYZ file reader 
  48.  
  49.     This reader takes a file in the XYZ Cartesian molecular model
  50.     format as specified by XMol (Copyright 1991 Research Equipment
  51.     Inc. dba Minnesota Supercomputer Center).  Three datasets,
  52.     which may be time-dependent, are then created: 
  53.     "name.atoms", "name.radii", and "name.form".  
  54.     The dataset, "name.atoms", may then be used to visualize/animate 
  55.     the molecule in ball-and/or-stick form using the Balls or Sticks
  56.     visualization objects.  The distance required for an inter-atomic bond
  57.     to be shown is determined by one of two methods:
  58.         1) the sum of the covalent radii, multiplied by a scale factor
  59.         2) an arbitrary (constant) value
  60.     As in the PDB reader, the size and colour of the balls and sticks are
  61.     determined by the atomic number & radius.  
  62.  
  63.    ---------------------------------- */
  64.  
  65. /* ----------------------------------
  66.     Converts XYZ name to atomic number
  67.     ---------------------------------- */
  68. #ifdef PROTO
  69. int AtomXYZNameToNumber(char *name)
  70. #else
  71. int AtomXYZNameToNumber(name)
  72. char *name;
  73. #endif
  74. /*Finds the atomic number for atom name.  Returns it or 0*/
  75. {
  76.  
  77.     int k;
  78.  
  79.  
  80.     for (k = 0; k < N_XYZ_SYMBOLS; ++k)
  81.     {
  82.  
  83.     
  84.     if (0 == strcmp(name, xyzSymbolTable[k] . name))
  85.     {
  86.     
  87.         return xyzSymbolTable[k] . atomicNumber;
  88.     }
  89.     }
  90.    /* if not found yet, check user-defined atoms */
  91.   for (k=0; k < numExtraXYZSymbols; k++)
  92.     {
  93.     if (0 == strcmp(name, xyzExtraSymbolTable[k] . name))
  94.     {
  95.  
  96.         return xyzExtraSymbolTable[k] . atomicNumber;
  97.     }
  98.     }
  99.  
  100.     return 0;
  101. }
  102.  
  103.  
  104. #define MISRAD            -1.0    /* Need to put one declaration somewhere*/
  105. #define STICKRADIUS            0.15
  106. #define MAXSTICKLENGTH         1.7
  107. #define CALC_COVALENTRADII     0    /* flag to calc bonds by covalent radii method */
  108. #define CALC_MAXSTICK        1    /* flag to calc bonds by arbitrary stick length */
  109. #define FUDGEFACTOR            1.0
  110.  
  111.  
  112. /*Ball style flags*/
  113. #define BS_JACKS    2
  114. #define BS_CUBES    4
  115. #define BS_OCTOHEDRA    8
  116. #define BS_SPHERES    16
  117.  
  118.  
  119.  
  120. /* -----------------------------------
  121.     XYZ file reader
  122.     modified from PDB file reader
  123.  
  124.    Format specification:
  125.  
  126.      "The XYZ format supports multi-step datasets.    Each  step  is
  127.       represented by a two-line "header," followed by one line for
  128.       each atom.
  129.  
  130.       The first line of a step's header is the number of atoms  in
  131.       that step.  This integer may be preceded by whitespace; any-
  132.       thing    on the line after the integer is ignored.  The    second
  133.       line    of  the     header     leaves    room for a descriptive string.
  134.       This line may    be blank, or it    may contain  some  information
  135.       pertinent to that particular step, but it _m_u_s_t exist,    and it
  136.       _m_u_s_t be just one line    long.
  137.  
  138.       Each line of text describing a single    atom must  contain  at
  139.       least     four  fields of information, separated    by whitespace:
  140.       the atom's type (a short string of alphanumeric characters),
  141.       and  its  x-,     y-, and z-positions. Optionally, extra    fields
  142.       may be used to specify a charge for the atom,    and/or a  vec-
  143.       tor  associated  with     the  atom.  If    an input line contains
  144.       five or eight    fields,    the fifth field    is interpreted as  the
  145.       atom's  charge;  otherwise, a    charge of zero is assumed.  If
  146.       an input line    contains seven or eight    fields,    the last three
  147.       fields are interpreted as the    components of a    vector.     These
  148.       components should be specified in angstroms."
  149.         (XYZ man pages, Minnesota Supercomputer Centre, 1991)
  150.     ----------------------------------- */
  151. #ifdef PROTO
  152. static ObjPtr ReadXYZFile(ObjPtr reader, char *name)
  153. #else
  154. static ObjPtr ReadXYZFile(reader,name)
  155. ObjPtr  reader;
  156. char *name;
  157. #endif
  158. {
  159.     /* -----------------------------------
  160.         Variables for dataset and data form
  161.         ----------------------------------- */
  162.        real curTime = 0.0;    /* keeps track of current animation frame */
  163.     char extension[5];        /* auto-numbering extension to dataset name*/
  164.                         /* if same dataset opened more than once */
  165.     /* NOTE: important to have new names for each dataset and form */
  166.     char atomsName[80], radiiName[80],  tmpName[80],formName[80];
  167.     char *tc;                /* terminating character */
  168.     real bounds[6];        /* xmin, xmax, ymin, .... */
  169.     long dims[2];            /* dimension vector for unstructured data form */
  170.     ObjPtr formVectors, dataForm, atomicNumberDataset, radiusDataset;
  171.     /* -----------------------------------
  172.         Variables for file input and parsing
  173.         ----------------------------------- */    
  174.         char cmdStr[256];
  175.         FILE *inFile;
  176.     real field5, field6, field7, field8;    /* store optional xyz fields of information*/
  177.     char s[4];
  178.     char atomNameBuf[4];
  179.     /* -----------------------------------
  180.         Variables for atoms
  181.         ----------------------------------- */    
  182.      long numAtoms;        /* number of atoms per frame */
  183.     atomData *currentAtoms;    /* points to array of data for current atoms read */
  184.     float centre[3];        /* coords of centre of atom */
  185.     real ballScaleFactor;
  186.     /* -----------------------------------
  187.         Variables for bonds
  188.         ----------------------------------- */    
  189.     long nBonds = 0;        /* number of bonds guessed */
  190.     TwoReals *bondsBuffer = 0;    
  191.     long nBondsAllocated = 0;
  192.     real maxStickLength;
  193.     ObjPtr var;
  194.     real  sum, temp;
  195.     real covalentRadius1, covalentRadius2, scaleFactor;
  196.     int bondCalcType;    /* method used to calculate minimum interatomic dist for bond*/
  197.  
  198.     long i, j,k, count;            /* counters*/
  199.     Bool convertToAngstroms;    /* false if already in angstroms, true otherwise */    
  200.  
  201.  
  202.     /* -------------------------------------
  203.       From the file reader controls, get:
  204.         -value of the arbitrary maximum bond length
  205.         -bond guessing method
  206.         -covalent radii fudge factor
  207.         -ball radius scale factor
  208.         -state of units conversion checkbox
  209.        ------------------------------------- */
  210.     var = GetRealVar("ReadXYZFile", reader, USER1);
  211.       if (!var)
  212.             maxStickLength = MAXSTICKLENGTH;
  213.       else
  214.             maxStickLength = GetReal(var);
  215.  
  216.     if ((var = GetIntVar("ReadXYZFile", reader, USER3)) == NULLOBJ)
  217.     {
  218.         bondCalcType = CALC_COVALENTRADII;
  219.     }
  220.     else
  221.     {
  222.         bondCalcType = GetInt(var);
  223.     }
  224.     if ((var = GetRealVar("ReadXYZFile", reader, USER4)) == NULLOBJ)
  225.     {
  226.         scaleFactor = FUDGEFACTOR;
  227.     }
  228.     else
  229.     {
  230.         scaleFactor = GetReal(var);
  231.     }
  232.      if ((var = GetRealVar("ReadXYZFile", reader, USER5)) == NULLOBJ)
  233.     {
  234.         ballScaleFactor = 1.0;
  235.     }
  236.     else
  237.     {
  238.         ballScaleFactor = GetReal(var);
  239.     }     
  240.     convertToAngstroms = GetPredicate(reader, USER6) ? true : false;
  241.  
  242.  
  243.     /* ------------------------------------ 
  244.       read the file in
  245.        ------------------------------------- */
  246.         inFile = fopen(name, "r");
  247.         if (!inFile)
  248.         {
  249.         Error("ReadXYZFile", OPENFILEERROR, name);
  250.         return NULLOBJ;
  251.         }
  252.  
  253.     /* -------------------------------------
  254.         make names for datasets 
  255.         ------------------------------------- */
  256.     strcpy (tmpName, name);
  257.     tc = tmpName;
  258.     while(*tc && (*tc != '@') && (*tc != '.'))
  259.         tc++;
  260.     *tc = '\0';
  261.  
  262.     strcpy(atomsName, tmpName);
  263.     strcpy(radiiName, tmpName);
  264.     strcpy(formName, tmpName);
  265.  
  266.     strcat(atomsName, ".atoms");
  267.     strcat(radiiName, ".radii");
  268.     strcat(formName, ".form");    
  269.  
  270.     k=1;
  271.     strcpy(tmpName, atomsName);
  272.     while (FindDatasetByName(tmpName))
  273.     {
  274.         ++k;
  275.         sprintf(tmpName, "%s(%d)", atomsName,k);
  276.     }
  277.         
  278.     if (k > 1)
  279.     {
  280.         sprintf(extension, "(%d)", k);
  281.         strcat(atomsName, extension);
  282.         strcat(radiiName, extension);
  283.         strcat(formName, extension);
  284.     }
  285.  
  286.  
  287.  
  288.     /* ----------------------------------------
  289.           make bounds impossible
  290.         ---------------------------------------- */
  291.     bounds[0] = bounds[2] = bounds[4] = PLUSINF;
  292.     bounds[1] = bounds[3] = bounds[5] = MINUSINF;
  293.  
  294.         /*------------------------------------------
  295.         Read the picture in
  296.            ------------------------------------------*/
  297.     while    (fgets(tempStr, 40, inFile)!= NULL)
  298.     {
  299.  
  300.         /* First line should have #atoms per frame */
  301.         if (1 == sscanf(tempStr, "%s", cmdStr))
  302.         {
  303.             if (1 != sscanf(tempStr, "%d", &numAtoms))
  304.             {
  305.                 FileFormatError("ReadXYZFile", "# atoms expected");
  306.                 return NULLOBJ;
  307.             }
  308.             currentAtoms = (atomData *) Alloc (numAtoms * sizeof(atomData));
  309.             if ( currentAtoms == NULL)
  310.             {
  311.                 FileFormatError("ReadXYZFile", "unable to allocate memory \
  312. for atoms");
  313.                 return NULLOBJ;
  314.             }
  315.         }
  316.         else
  317.         {
  318.             FileFormatError("ReadXYZFile", "# atoms expected");
  319.             return NULLOBJ;
  320.         }
  321.  
  322.         /* Skip the second line, which should be just a comment */
  323.         ReadLn(inFile);
  324.         
  325.         /* Read a line in for each atom */
  326.         for (i = 0; i < numAtoms; i++ )
  327.         {
  328.             fgets( tempStr, 80, inFile);
  329.  
  330.             /* ----------------------------------
  331.                 format of each line may be one of:
  332.  
  333.                 "name x y z\n"
  334.                 "name x y z charge\n"
  335.                 "name x y z vectorX vectorY vectorZ\n"
  336.                 "name x y z charge vectorX vectorY vectorZ\n"
  337.                ---------------------------------- */
  338.             if (sscanf(tempStr, "%s %g %g %g %g %g %g %g", &s, &(centre[0]), \
  339.             &(centre[1]),  &(centre[2]), &field5, &field6, &field7, &field8) \
  340.             >= 4 )
  341.             {
  342.                 /* -----------------------------
  343.                     Record atomic name and number
  344.                     and convert units if necessary
  345.                    ------------------------------ */
  346.                 strcpy(currentAtoms[i].name, s);
  347.                 currentAtoms[i].atomicNumber = AtomXYZNameToNumber(currentAtoms[i].name);
  348.                 
  349.                     
  350.  
  351.                 for (j = 0; j < 3; j++ )
  352.                 {
  353.                     if (convertToAngstroms)
  354.                         centre[j] *= ANGSTROMS_PER_ALU;
  355.  
  356.                     /* store all atom coords for later */
  357.                     currentAtoms[i].coords[j] = centre[j];
  358.                     
  359.                     /* revise bounds if necessary */
  360.                     if ( centre[j] < bounds[j*2] )
  361.                         bounds[j*2] = centre[j];
  362.                     if ( centre[j] > bounds[j*2 + 1] )
  363.                         bounds[j*2 + 1] = centre[j];
  364.             
  365.                 }        
  366.  
  367.             }
  368.             else
  369.             {
  370.                 FileFormatError("ReadXYZFile", "Badly formatted atom");
  371.             }
  372.         }
  373.         nBonds = 0;
  374.         nBondsAllocated = 200;
  375.         bondsBuffer = (TwoReals *) Alloc(sizeof(TwoReals) * nBondsAllocated);    
  376.         switch ( bondCalcType )
  377.         {
  378.             case CALC_COVALENTRADII:
  379.                 for ( i = 0; i < (numAtoms-1); i++ )
  380.                 {
  381.                     covalentRadius1 = atomInfo[currentAtoms[i].atomicNumber-1].covalentRadius;    
  382.                     /* --------------------------------
  383.                         If either atom is missing the covalent radius data,
  384.                         do not draw a bond between them.
  385.                         -------------------------------- */
  386.  
  387.                     if (covalentRadius1 != MISRAD)
  388.                     {
  389.                         for (j = (i+1); j < numAtoms; j++ )
  390.                         {
  391.                             covalentRadius2 = atomInfo[currentAtoms[j].atomicNumber-1].covalentRadius;
  392.  
  393.                             if ( covalentRadius2 != MISRAD)
  394.                             {
  395.                                 sum = 0.0;
  396.                                 for (k = 0; k < 3; k++)
  397.                                 {
  398.                                     temp = currentAtoms[j].coords[k] - \
  399.                                         currentAtoms[i].coords[k];
  400.                                     sum += temp * temp;
  401.                                 }
  402.                                 if (sum <= SQUARE ( (covalentRadius1 + covalentRadius2)*scaleFactor))
  403.                                 {
  404.                                         if (nBonds >= nBondsAllocated)
  405.                                         {
  406.                                         nBondsAllocated += 200;
  407.                                         bondsBuffer = (TwoReals *) Realloc(bondsBuffer,
  408.                                         sizeof(TwoReals) * nBondsAllocated);
  409.                                         }
  410.                                         bondsBuffer[nBonds][0] = i;
  411.                                         bondsBuffer[nBonds][1] = j;
  412.                                         ++nBonds;
  413.                                 }    
  414.                             }
  415.                         }            
  416.                     }
  417.                 }
  418.                 break;
  419.             case CALC_MAXSTICK:
  420.                 for ( i = 0; i < (numAtoms-1); i++ )
  421.                 {
  422.                     for (j = (i+1); j < numAtoms; j++ )
  423.                     {
  424.                         sum = 0.0;
  425.                         for (k = 0; k < 3; k++)
  426.                         {
  427.                             temp = currentAtoms[j].coords[k] - \
  428.                                     currentAtoms[i].coords[k];
  429.                             sum += temp * temp;
  430.                         }
  431.     
  432.                         if (sum <= SQUARE(maxStickLength))
  433.                         {
  434.                                 if (nBonds >= nBondsAllocated)
  435.                                 {
  436.                                 nBondsAllocated += 200;
  437.                                 bondsBuffer = (TwoReals *) Realloc(bondsBuffer,
  438.                                 sizeof(TwoReals) * nBondsAllocated);
  439.                                 }
  440.                                 bondsBuffer[nBonds][0] = i;
  441.                                 bondsBuffer[nBonds][1] = j;
  442.                                 ++nBonds;
  443.                         }
  444.                     }
  445.                 }
  446.                 break;
  447.         }
  448.  
  449.  
  450.         /* make datasets */
  451.         atomicNumberDataset = NewDataset( atomsName, 1, &numAtoms, 0);
  452.  
  453.         /* this variable signals the palette creation later */
  454.         SetVar(atomicNumberDataset, UNITSNAME, NewString("Atomic number"));
  455.         
  456.         radiusDataset = NewDataset( radiiName, 1, &numAtoms, 0 );
  457.  
  458.         /* --------------------------------------
  459.             make a vector dataset for the form
  460.             -------------------------------------- */
  461.         formVectors = NewDataset(formName, 1, &numAtoms, 3);
  462.         dims[0] = numAtoms;
  463.         dims[1] = nBonds;
  464.         if (nBonds > 0)
  465.         { /* one dimensional dataset (points and lines)  */
  466.             dataForm = NewUnstructuredDataForm(tmpName, 1, dims, bounds,
  467. formVectors);
  468.         }
  469.         else    
  470.         {/* zero dimensional dataset (points only)*/
  471.             dataForm = NewUnstructuredDataForm(tmpName, 0, dims, bounds,
  472. formVectors);
  473.         }
  474.  
  475.         SetDatasetForm(atomicNumberDataset, dataForm);
  476.         SetDatasetForm(radiusDataset, dataForm);
  477.         /* Fill it all up */
  478.         if (nBonds)
  479.         {
  480.         /*Set the bonds*/
  481.         for (k = 0; k < nBonds; ++k)
  482.         {
  483.             Put1Edge(dataForm, k,
  484.             (long) bondsBuffer[k][0], (long) bondsBuffer[k][1]);
  485.         }
  486.         }
  487.         SetCurField(FORMFIELD, formVectors);
  488.         SetCurField(FIELD1, atomicNumberDataset);
  489.         SetCurField(FIELD2, radiusDataset);
  490.         
  491.         for (i=0; i < numAtoms; i++ )
  492.         {
  493.                 PutFieldComponent(FORMFIELD, 0, &i, currentAtoms[i].coords[0]);
  494.                 PutFieldComponent(FORMFIELD, 1, &i, currentAtoms[i].coords[1]);
  495.                 PutFieldComponent(FORMFIELD, 2, &i, currentAtoms[i].coords[2]);
  496.     
  497.             
  498.  
  499.                 PutFieldComponent(FIELD1, 0, &i, currentAtoms[i].atomicNumber ? (real)
  500. currentAtoms[i].atomicNumber : missingData);
  501.                 PutFieldComponent(FIELD2, 0, &i, currentAtoms[i].atomicNumber ?
  502. atomInfo[currentAtoms[i].atomicNumber-1].radius : missingData );
  503.         }
  504.  
  505.         /* -------------------------
  506.             create new variables of SIZEOBJ ID
  507.             allows sizing and scaling of atomic radii
  508.            -------------------------- */
  509.         SetVar(atomicNumberDataset, SIZEOBJ, radiusDataset);
  510.         SetVar(radiusDataset, SIZEOBJ, radiusDataset);
  511.  
  512.         SetVar(atomicNumberDataset, POINTSTYLE, NewInt(BS_SPHERES));
  513.  
  514.         /* ------------------------------
  515.             Interpolation in time
  516.              ------------------------------ */
  517.  
  518.         SetVar(atomicNumberDataset, INTERPOLATEP, ObjTrue);
  519.         SetVar(radiusDataset, INTERPOLATEP, ObjTrue);
  520.         SetVar(formVectors, INTERPOLATEP, ObjTrue);
  521.  
  522.         RegisterTimeDataset ( atomicNumberDataset, curTime );
  523.         RegisterTimeDataset ( radiusDataset, curTime );
  524.         /* yes I know radii don't change with time */
  525.         RegisterTimeDataset ( formVectors, curTime );
  526.  
  527.          curTime += 1.0;    /* increment frame count */
  528.         Free(currentAtoms);
  529.     
  530.     }
  531.  
  532.  
  533.     fclose(inFile);
  534.   
  535.     /* ---------------------------------
  536.         HACK!
  537.  
  538.         Sets default ball visualization to a sphere instead
  539.         of an octohedron.  SetVar of POINTSTYLE of 
  540.         dataset itself doesn't seem to work.
  541.     ----------------------------------- */
  542.       SetVar(visBalls, POINTSTYLE, NewInt(BS_SPHERES));
  543.     /* ---------------------------------
  544.         MORE HACK!
  545.  
  546.         Set scale factor of balls' size object
  547.         ---------------------------------- */
  548.     SetVar(visBalls, SIZEFACTOR, NewReal(ballScaleFactor));
  549.  
  550.         return NULLOBJ;
  551. }
  552.  
  553. /* ------------------------------------
  554.     Read in extra user-defined symbolic information for XYZ file format
  555.     ------------------------------------ */
  556. #ifdef PROTO
  557. int     ReadExtraXYZSymbolTable (void)
  558. #else
  559. int     ReadExtraXYZSymbolTable()
  560. #endif
  561. {
  562.     char tempStr[80];
  563.     int i, numTypes;
  564.     FILE *inFile;
  565.     char symbol[3], name[6];
  566.     int number;
  567.     
  568.     /* ------------------------------------ 
  569.       read the file in
  570.        ------------------------------------- */
  571.         inFile = fopen(xyzExtraInfoFile, "r");
  572.         if (!inFile)
  573.         {
  574.         /* no extra symbols */        
  575.         return(0);
  576.         }
  577.  
  578.     /* -------------------------------------
  579.         first line should have number of types listed in file 
  580.         ------------------------------------- */
  581.     if (fgets(tempStr, 80, inFile) == NULL)
  582.         return(-1);
  583.         
  584.     if (sscanf(tempStr, "%d", &numTypes) != 1)
  585.         return(-1);
  586.     
  587.     /* --------------------------------------
  588.         allocate memory to store symbol table
  589.         -------------------------------------- */
  590.     xyzExtraSymbolTable = (xyzInfo *) Alloc ( numTypes * sizeof (xyzInfo) );
  591.     if (xyzExtraSymbolTable == NULL)
  592.         return (-1);
  593.  
  594.     for ( i = 0; i < numTypes; i++ )
  595.     {
  596.         if (fgets(tempStr, 80, inFile) == NULL)
  597.             return (-1);
  598.  
  599.         if (sscanf (tempStr, "%s %s %d", symbol, name, &number ) == 3 )
  600.         {
  601.             strcpy(xyzExtraSymbolTable[i].symbol,symbol);
  602.             strcpy(xyzExtraSymbolTable[i].name, name);
  603.             xyzExtraSymbolTable[i].atomicNumber = number;
  604.         }
  605.         else
  606.             return(-1);
  607.     }
  608.  
  609.     return (numTypes);    
  610. }
  611.  
  612.  
  613. /* --------------------------------
  614.     XYZ File Reader Controls 
  615.     Modified from AddNXRControls()
  616.    -------------------------------- */
  617. #define XYZ_ITEMWIDTH 225
  618. #define XYZ_MIDDLE_SEP 45
  619. #define XYZ_SMALLER_ITEM 180
  620. #define XYZ_RADIOWIDTH 400
  621. #define XYZ_BOXLEFT 325
  622. #define XYZ_BOXWIDTH 60
  623.  
  624. #ifdef PROTO
  625. static ObjPtr    AddXYZControls(ObjPtr FileReader, ObjPtr PanelContents)
  626. #else
  627. static ObjPtr    AddXYZControls(FileReader, PanelContents)
  628.   ObjPtr  FileReader, PanelContents;
  629. #endif
  630. {
  631.   ObjPtr titleBox, radio, checkBox, button;
  632.   ObjPtr  checkbox, textbox, var;
  633.   int     left, right, bottom, top;
  634.  
  635.   bottom = left = MAJORBORDER;
  636.  
  637.     /*Create the check box for units conversion*/
  638.     top = bottom + CHECKBOXHEIGHT;
  639.     right = left + FRTIMEWIDTH;
  640.     checkBox = NewCheckBox(left, right, bottom, top,
  641.         "Coordinates given in atomic length units",
  642.         GetPredicate(FileReader, USER6));
  643.     SetVar(checkBox, PARENT, PanelContents);
  644.     PrefixList(PanelContents, checkBox);
  645.     bottom = top + MINORBORDER;
  646.     if (!GetVar(FileReader, USER6))
  647.     {
  648.         SetVar(FileReader, USER6, ObjFalse);
  649.     }
  650.     AssocDirectControlWithVar(checkBox, FileReader, USER6);
  651.     SetVar(checkBox, HELPSTRING, 
  652.         NewString("If this box is checked, the file reader will consider all \
  653. xyz coordinates to be in ALU and convert all coordinates to Angstroms. \
  654. A conversion factor of 0.5292 Angstroms per ALU is used in the \
  655. calculation."));
  656.  
  657.  
  658.  /* Pre-set balls scaling factor textbox*/
  659.   left = MAJORBORDER;
  660.   right = left + XYZ_SMALLER_ITEM;
  661.   bottom = top;
  662.   top = bottom + EDITBOXHEIGHT;
  663.  
  664.   textbox = NewTextBox(left, right, bottom, top,
  665.                0, "Ball Radius Scaling Factor Text", "Ball Radius Scaling Factor:");
  666.   SetVar(textbox, PARENT, PanelContents);
  667.   PrefixList(PanelContents, textbox);
  668.  
  669.     left = right;
  670.   right = left + XYZ_BOXWIDTH;
  671.  
  672.  
  673.   var = GetRealVar("AddXYZControls", FileReader, USER5);
  674.   if (!var)
  675.     SetVar(FileReader, USER5, NewReal(1.0));
  676.   textbox = NewTextBox(left, right, bottom, top,
  677.                EDITABLE + WITH_PIT + ONE_LINE,
  678.                "Ball Radius Scaling Factor", "0");
  679.   SetVar(textbox, PARENT, PanelContents);
  680.   PrefixList(PanelContents, textbox);
  681.   SetTextAlign(textbox, RIGHTALIGN);
  682.   AssocTextRealControlWithVar(textbox, FileReader, USER5, 
  683.                   0, plusInf, TR_NE_TOP);
  684.   SetVar(textbox, HELPSTRING,
  685.      NewString("This value indicates the scaling factor applied to the radius of \
  686. each atom (ball)  in ball-and-stick molecular visualization.  This option allows \
  687. the user to pre-set the scale factor of a size object in a balls visualization.  For \
  688. example, a user reading in a lot of large organic molecules may want to set \
  689. this value to, say, 0.2."));
  690.  
  691.  
  692.  
  693.  
  694.     left = MAJORBORDER;
  695.     bottom = top + MINORBORDER;
  696.       /*Do the title box around the radio group*/
  697.     right = left + 2 * MINORBORDER + XYZ_RADIOWIDTH;
  698.     top = bottom + 2 * MINORBORDER + 1 * CHECKBOXSPACING + 2 \
  699. * CHECKBOXHEIGHT + TITLEBOXTOP;
  700.  
  701.     titleBox = NewTitleBox(left, right, bottom, top, "Calculate Molecular Bonds");
  702.     PrefixList(PanelContents, titleBox);
  703.  
  704.     SetVar(titleBox, PARENT, PanelContents);
  705.  
  706.     /*Make the radio buttons*/
  707.     radio = NewRadioButtonGroup("Molecular Bonds Radio");
  708.     SetVar(radio, PARENT, PanelContents);
  709.     SetVar(radio, REPOBJ, FileReader);
  710.     PrefixList(PanelContents, radio);
  711.  
  712.     left += MINORBORDER;
  713.     right -= MINORBORDER;
  714.     top -= TITLEBOXTOP + MINORBORDER;
  715.     bottom = top - CHECKBOXHEIGHT;
  716.     button = NewRadioButton(left, right, bottom, top, "Sum of covalent radii * scale factor of:");
  717.     AddRadioButton(radio, button);
  718.     SetVar(button, HELPSTRING, NewString("When this button is down, the minimum \
  719. interatomic distance required for a bond (stick) to be drawn is calculated as the sum of the \
  720. respective covalent radii as defined in ScianPeriodicTable.h multiplied by a scale factor \
  721. (default 1.0)."));
  722.  
  723.     /* scale factor text box */
  724.   left = MAJORBORDER + XYZ_BOXLEFT;
  725.   right = left + XYZ_BOXWIDTH;
  726.   top = bottom + EDITBOXHEIGHT;
  727.  
  728.  
  729.   var = GetRealVar("AddXYZControls", FileReader, USER4);
  730.   if (!var)
  731.     SetVar(FileReader, USER4, NewReal(FUDGEFACTOR));
  732.   textbox = NewTextBox(left, right, bottom, top,
  733.                EDITABLE + WITH_PIT + ONE_LINE,
  734.                "Scale factor", "0");
  735.   SetVar(textbox, PARENT, PanelContents);
  736.   PrefixList(PanelContents, textbox);
  737.   SetTextAlign(textbox, RIGHTALIGN);
  738.  
  739.   AssocTextRealControlWithVar(textbox, FileReader, USER4, 
  740.                   0, plusInf, TR_NE_TOP);
  741.   SetVar(textbox, HELPSTRING,
  742.      NewString("This value is the scale factor by which the sum of covalent radii is multiplied \
  743. to determine minimum separation distance required between atoms for a bond (stick) to be drawn. \
  744. I.e. if this value is SF, then the atoms must be at most (covalent_radius1 + covalent_radius2)*SF \
  745. apart for a bond to occur."));
  746.  
  747.  
  748.     top = bottom - CHECKBOXSPACING;
  749.     bottom = top - EDITBOXHEIGHT;
  750.     top = bottom + CHECKBOXHEIGHT;
  751.     left = MAJORBORDER+MINORBORDER;
  752.     button = NewRadioButton(left, right, bottom, top, "Arbitrary maximum stick length of:");
  753.     AddRadioButton(radio, button);
  754.     SetVar(button, HELPSTRING, NewString("When this button is down, the minimum \
  755. interatomic distance required for a bond (stick) to be drawn is simply the maximum \
  756. stick length."));
  757.  
  758.   /* Maximum stick length setting text box*/
  759.   left = MAJORBORDER + XYZ_BOXLEFT;
  760.   right = left + XYZ_BOXWIDTH;
  761.   top = bottom + EDITBOXHEIGHT;
  762.  
  763.  
  764.   var = GetRealVar("AddXYZControls", FileReader, USER1);
  765.   if (!var)
  766.     SetVar(FileReader, USER1, NewReal(MAXSTICKLENGTH));
  767.   textbox = NewTextBox(left, right, bottom, top,
  768.                EDITABLE + WITH_PIT + ONE_LINE,
  769.                "Maximum stick length", "0");
  770.   SetVar(textbox, PARENT, PanelContents);
  771.   PrefixList(PanelContents, textbox);
  772.   SetTextAlign(textbox, RIGHTALIGN);
  773.  
  774.   AssocTextRealControlWithVar(textbox, FileReader, USER1, 
  775.                   0, plusInf, TR_NE_TOP);
  776.   SetVar(textbox, HELPSTRING,
  777.      NewString("This value indicates the minimum separation distance required between \
  778. atoms for a bond (stick) to be drawn."));
  779.  
  780.     
  781.     var = GetIntVar("AddXYZControls", FileReader, USER3);
  782.     if (!var) SetVar(FileReader, USER3, NewInt(CALC_COVALENTRADII));
  783.     AssocDirectControlWithVar(radio, FileReader, USER3);
  784.     SetVar(radio, HELPSTRING, NewString("This radio button group controls how \
  785. the minimum interatomic distance required for a bond is calculated."));
  786.   
  787.  
  788.  
  789.  
  790.   return NULLOBJ;
  791. }
  792.  
  793.  
  794. /* -------------------------------------------
  795.     GRD file reader 
  796.  
  797.     This file reader takes the 3D grid data produced by the PLOT
  798.     keyword of BioSym's software, DMol, as described on pages
  799.     B-31 to B-35 of the DMol User Manual.
  800.     ------------------------------------------- */
  801. #ifdef PROTO
  802. static int parseGRDHeader(FILE *infile, gridParameters *p)
  803. #else
  804. static int parseGRDHeader(infile, p)
  805. FILE *infile;
  806. gridParameters *p;
  807. #endif
  808. {
  809.     char tempStr[80];
  810.     int index;
  811.  
  812.     /* -----------------------------
  813.         Parse header, expecting exactly 5 lines 
  814.         ----------------------------- */
  815.     ReadLn(infile);        /* first line is a comment; ignore */
  816.  
  817.     ReadLn(infile);        /* second line is a FORTRAN format description
  818.                             should be 1F15.10 (one data item per line)
  819.                             ignore for now */
  820.  
  821.     if (fgets(tempStr, 80, infile) != NULL)    /* third line has cell lengths and angles*/
  822.     {
  823.         if (6 != sscanf( tempStr, "%g %g %g %g %g %g", &p->a,&p->b,&p->c, \
  824.                                   &p->alpha, &p->beta, &p->gamma) )
  825.         {
  826.             FileFormatError("ReadGRDFile", "Bad format for cell lengths and angles");
  827.             return(-1);
  828.         }
  829.     }
  830.     else
  831.     {
  832.         FileFormatError("ReadGRDFile", "Expected cell lengths and angles");
  833.         return(-1);
  834.     }
  835.     if (fgets(tempStr, 80, infile) != NULL)    /* fourth line has #grid spaces*/
  836.     {
  837.         if (3 != sscanf(tempStr, "%d %d %d", &p->numX, &p->numY, &p->numZ))
  838.         {
  839.             FileFormatError("ReadGRDFile", "Bad format for grid spaces");
  840.             return(-1);
  841.         }
  842.     }
  843.     else
  844.     {
  845.         FileFormatError("ReadGRDFile", "Expected grid space information");
  846.         return(-1);
  847.     }
  848.  
  849.     if (fgets(tempStr, 80, infile) != NULL)    /* fifth line has fastest varying index*/
  850.     {                                    /* and grid spaces to either side of origin*/
  851.         if (7 != sscanf(tempStr, "%d %d %d %d %d %d %d", &(index), &p->xStart, \
  852.         &p->xEnd, &p->yStart, &p->yEnd, &p->zStart, &p->zEnd) )
  853.         {
  854.             FileFormatError("ReadGRDFile", "Bad format for fastest varying index or \
  855. grid spaces to either side of origin");
  856.             return(-1);
  857.         }
  858.         if ( index == 1)
  859.             p->xIsFastest = true;
  860.         else if (index == 3)
  861.             p->xIsFastest = false;
  862.         else
  863.         {
  864.             FileFormatError("ReadGRDFile", "Fastest varying index should be either \
  865. 1 or 3");
  866.             return(-1);
  867.         }
  868.     }
  869.     else
  870.     {
  871.         FileFormatError("ReadGRDFile", "Expected fastest varying index and grid \
  872. spaces to either side of origin information");
  873.         return(-1);
  874.     }
  875. }
  876.  
  877. #ifdef PROTO
  878. static ObjPtr ReadGRDFile(ObjPtr reader, char *name)
  879. #else
  880. static ObjPtr ReadGRDFile(reader,name)
  881. ObjPtr  reader;
  882. char *name;
  883. #endif
  884. /*Reads a GRD file from name*/
  885. {
  886.     FILE *infile;
  887.         ObjPtr retVal;                
  888.     ObjPtr dataForm;
  889.     ObjPtr gridDataset;
  890.       long    dims[3];
  891.       real    bounds[6];    /* Used to define parameters of datasets */
  892.       long    indices[3];    /* and dataforms */
  893.     char formName[80], dataName[80];
  894.     char *tc;
  895.     gridParameters grid;
  896.     real curData;
  897.     char tempStr[40];
  898.     real min, max, minMax[2];
  899.     ObjPtr minMaxArray;    
  900.  
  901.       /*---------------------
  902.          Derive a dataset name from the file name 
  903.        --------------------- */
  904.  
  905.       strcpy(formName, name);
  906.       tc = formName;
  907.       while (*tc && (*tc != '@') && (*tc != '.'))
  908.         tc++;
  909.       *tc = '\0';
  910.       strcpy(dataName, formName);
  911.  
  912.       strcat(formName, ".form");
  913.       strcat(dataName, ".data");
  914.   
  915.  
  916.   /* Open the file and verify that it contains appropriate data */
  917.  
  918.       if ((infile = fopen(name, "r")) == NULL) 
  919.     {
  920.             FileFormatError("ReadGRDFile", "File not accessible");
  921.             return NULLOBJ;
  922.       }
  923.  
  924.     if (parseGRDHeader(infile, &grid) == -1)
  925.     {
  926.         FileFormatError("ReadGRDFile", "Unable to parse grid parameters");
  927.         return NULLOBJ;
  928.     }
  929.  
  930.     /* -----------------------------------
  931.         Calculate the dimensions of the grid.
  932.         There is always one more grid point than number
  933.         of grid spaces in each direction.
  934.        ------------------------------------ */
  935.     dims[0] = grid.numX + 1;
  936.     dims[1] = grid.numY + 1;
  937.     dims[2] = grid.numZ + 1;
  938.  
  939.     /* -----------------------------------
  940.         Create the scalar dataset
  941.         ----------------------------------- */
  942.     retVal = NewStructuredDataset(dataName, 3, dims, 0 );
  943.  
  944.     /* ------------------------------------
  945.         Calculate grid bounds 
  946.         ------------------------------------ */
  947.     bounds[0] =  grid.a * (real)grid.xStart / ((real)grid.xEnd-(real)grid.xStart); /*xMin */
  948.     bounds[1] = grid.a * ((real) (grid.xStart + grid.numX))/((real) (grid.xEnd-grid.xStart));
  949.     bounds[2] =  grid.b * (real)grid.yStart / ((real)grid.yEnd-(real)grid.yStart); /*yMin */
  950.     bounds[3] = grid.b * ((real) (grid.yStart + grid.numY))/((real) (grid.yEnd-grid.yStart));
  951.     bounds[4] =  grid.c * (real)grid.zStart / ((real)grid.zEnd-(real)grid.zStart); /*zMin */
  952.     bounds[5] = grid.c * ((real) (grid.zStart + grid.numZ))/((real) (grid.zEnd-grid.zStart));
  953.  
  954.     /* ------------------------------------
  955.         Create a regular data form
  956.         ------------------------------------  */
  957.     dataForm = NewRegularDataForm(formName, 3, dims, bounds);
  958.     SetDatasetForm(retVal, dataForm);
  959.     /* ------------------------------------
  960.         Fill up the dataset
  961.         depending on which index varies fastest
  962.         Keep track of minimum and maximum values as
  963.         data values are read in.
  964.         ------------------------------------ */
  965.     SetCurField(FIELD1, retVal);  /* Link the field to the data form */
  966.     min = PLUSINF;                /* Initialize min and max to */
  967.     max = MINUSINF;           /* extreme values */    
  968.  
  969.     if (grid.xIsFastest == true)
  970.     {
  971.         for (indices[2] = 0; indices[2] < dims[2]; (indices[2])++)
  972.         {
  973.             for (indices[1] = 0; indices[1] < dims[1]; (indices[1])++)
  974.             {
  975.                 for (indices[0] = 0; indices[0] < dims[0]; (indices[0])++)
  976.                 {
  977.                     if (fgets(tempStr, 40, infile) != NULL)
  978.                     {
  979.                         if (sscanf(tempStr, "%g", &curData) == 1)
  980.                         {
  981.                             PutFieldComponent(FIELD1, 0, indices,curData);
  982.                             if (curData < min)
  983.                                 min = curData;
  984.                             if (curData > max)
  985.                                 max = curData;
  986.  
  987.                         }
  988.                         else
  989.                         {
  990.                             FileFormatError("ReadGRDFile", "poor data format");
  991.                             return(NULLOBJ);
  992.                         }
  993.                     }
  994.                     else
  995.                     {
  996.                         FileFormatError("ReadGRDFile", "expected more data");
  997.                         return(NULLOBJ);
  998.                     }
  999.                 }
  1000.             }
  1001.         }
  1002.     }
  1003.     else    /* y varies fastest */
  1004.     {
  1005.         for (indices[2] = 0; indices[2] < dims[2]; (indices[2])++)
  1006.         {
  1007.             for (indices[0] = 0; indices[0] < dims[0]; (indices[0])++)
  1008.             {
  1009.                 for (indices[1] = 0; indices[1] < dims[1]; (indices[1])++)
  1010.                 {
  1011.                     if (fgets(tempStr, 40, infile) != NULL)
  1012.                     {
  1013.                         if (sscanf(tempStr, "%g", &curData) == 1)
  1014.                         {
  1015.                             PutFieldComponent(FIELD1, 0, indices,curData);
  1016.  
  1017.                             if (curData < min)
  1018.                                 min = curData;
  1019.                             if (curData > max)
  1020.                                 max = curData;
  1021.  
  1022.                         }
  1023.                         else
  1024.                         {
  1025.                             FileFormatError("ReadGRDFile", "poor data format");
  1026.                             return(NULLOBJ);
  1027.                         }
  1028.                     }
  1029.                     else
  1030.                     {
  1031.                         FileFormatError("ReadGRDFile", "expected more data");
  1032.                         return(NULLOBJ);
  1033.                     }
  1034.                 }
  1035.             }
  1036.         }
  1037.     }
  1038.     /* -------------------------
  1039.         Record the minimum and maximum data values
  1040.         for colour palette and isosurface scales 
  1041.         ------------------------- */
  1042.  
  1043.     minMax[0] = min;
  1044.     minMax[1] = max; 
  1045.     minMaxArray = NewRealArray(1, 2L);
  1046.     CArray2Array(minMaxArray, minMax);
  1047.     SetVar(retVal, MINMAX, minMaxArray);
  1048.  
  1049.     RegisterDataset(retVal);
  1050.     fclose(infile);
  1051.     return(NULLOBJ);
  1052. }
  1053.  
  1054. /* ------------------------------------
  1055.     Gaussian File Readers: GSN1, GSN2, GSN
  1056.  
  1057.     The purpose of each file reader is to take as input
  1058.     the coefficients and exponents of gaussian primitives,
  1059.     obtained from ab initio calculations such as that produced
  1060.     by G90, and approximate the molecular orbitals by a linear
  1061.     combination of the primitives. 
  1062.  
  1063.     The value computed at each grid point is obtained from
  1064.     the function:
  1065.  
  1066.     f(x,y,z) = SUM(i=1..n){ d_i * c_i * g_i (zeta,x,y,z,nx,ny,nz,
  1067.                                               Rx, Ry, Rz ) }
  1068.     where:
  1069.         d_i:    molecular orbital coefficient
  1070.         c_i: S,P,D shell coefficient of basis function
  1071.         g_i: unnormalized Cartesian Gaussian function
  1072.             with origin at R
  1073.         zeta: orbital exponent
  1074.         nx,ny,nz: determine type of shell
  1075.   
  1076.     (see "Simplified Introduction to Ab Initio Basis Sets.
  1077.     Terms and Notation", by Jan K. Labanowski, Ohio
  1078.     SuperComputer Center.  email: jkl@osc.edu,
  1079.     JKL@OHSTPY.BITNET
  1080.     see also Obara and Saika, "Efficient recursive computation 
  1081.     of molecular integrals over Cartesian Gaussian functions",
  1082.     J.Chem. Phys. 84(7) (April 1986) p. 3964 )
  1083.  
  1084.     GSN1 and GSN2 format are essentially the same format
  1085.     (8-column and 9-column, respectively) for output from
  1086.     programs other than G90.  The files specify grid dimensions
  1087.     and the parameters of each term in the summation.
  1088.  
  1089.     GSN format is for G90 output.  The parameters are 
  1090.     extracted from G90 output for a particular orbital number,
  1091.     the molecular orbital computed, and the molecule 
  1092.     visualized in ball-and-stick form.
  1093.     ------------------------------------ */        
  1094.  
  1095. /* ------------------------------------
  1096.     Cartesian Gaussian function
  1097.     ------------------------------------ */
  1098. #ifdef PROTO
  1099. static real gaussianFunction(real x, real y, real z, real Rx, real Ry, real Rz, \
  1100. real zeta, int nx, int ny, int nz, real N)
  1101. #else
  1102. static real gaussianFunction( x, y, z, Rx, Ry, Rz, zeta, nx, ny, nz, N )
  1103. real x, y, z, Rx, Ry, Rz, zeta;
  1104. int nx, ny, nz;
  1105. real N;    /* normalization constant */
  1106. #endif
  1107. {    
  1108.     double term1, term2, term3, term4;
  1109.     double dist;
  1110.     real ans;
  1111.     int i;
  1112.  
  1113.     term1 = 1.0;
  1114.     for (i=0; i < nx; i++ )
  1115.         term1 *= (x-Rx);
  1116.  
  1117.     term2 = 1.0;
  1118.     for (i=0; i < ny; i++)
  1119.         term2 *= y-Ry;
  1120.     
  1121.     term3 = 1.0;
  1122.     for (i=0; i < nz; i++)
  1123.         term3 *= z-Rz;
  1124.  
  1125.     dist = (double) ( (x-Rx)*(x-Rx) + (y-Ry)*(y-Ry) + (z-Rz)*(z-Rz) );
  1126.     term4 = exp ( (double) ( -zeta * dist ) );
  1127.     ans = (real) ( ((double) N) * term1 * term2 * term3 * term4 );
  1128.     return(ans);
  1129. }
  1130. /* -----------------------
  1131.     oddSkipFactorial(n) = n!!
  1132.  
  1133.     returns 1 * 3 * 5 * 7 * .... n for n > 0, n odd
  1134.     returns 1 for n = -1
  1135.     ----------------------- */
  1136. #ifdef PROTO
  1137. static int oddSkipFactorial( int n )
  1138. #else
  1139. static int oddSkipFactorial(n)
  1140. int n;
  1141. #endif
  1142. {
  1143.     int k, ans;
  1144.  
  1145.     /* assume that n is odd */
  1146.     if (n < 0 )
  1147.     {
  1148.         if (n == -1)
  1149.             return(1);
  1150.         else
  1151.             return(0);
  1152.     }
  1153.  
  1154.     ans = 1;
  1155.     for ( k = 1; k <= n; k += 2)
  1156.         ans *= k;
  1157.     return(ans);
  1158. }
  1159.  
  1160. #ifdef PROTO
  1161. static real normalizationConstant ( real zeta, int nx, int ny, int nz )
  1162. #else
  1163. static real normalizationConstant ( zeta, nx, ny, nz )
  1164. real zeta;
  1165. int nx, ny, nz;
  1166. #endif
  1167. {
  1168.     double term1, term2, term3;
  1169.     double powTerm2, fact1, fact2, fact3;
  1170.     double inversePi = M_1_PI;
  1171.     real ans;
  1172.  
  1173.     term1 = pow( 2.0 * (double) zeta * inversePi , 0.75 );
  1174.     powTerm2 = ((double) (nx + ny + nz)) * 0.5;
  1175.     term2 = pow( 4.0 * (double) zeta, powTerm2 );
  1176.     fact1 = oddSkipFactorial( 2*nx - 1 );
  1177.     fact2 = oddSkipFactorial(2*ny - 1);
  1178.     fact3 = oddSkipFactorial(2*nz - 1);
  1179.     term3 = pow( (double) (fact1 * fact2 * fact3), -0.5);
  1180.      
  1181.     ans = (real) (term1 * term2 * term3 );
  1182.     return(ans);
  1183. }
  1184.  
  1185. /* --------------------------------------
  1186.     ParseGSN12Header
  1187.     
  1188.     Gets grid bounds, dimensions, and number of terms
  1189.     in summation for GSN1 and GSN2 file readers
  1190.    -------------------------------------- */
  1191. #ifdef PROTO
  1192. static int ParseGSN12Header( FILE *infile, real *px1, real *px2, real*py1, real *py2, real *pz1, \
  1193. real *pz2, int *pnumx, int *pnumy, int *pnumz, int *pn )
  1194. #else
  1195. static int ParseGSN12Header ( infile, px1, px2, py1, py2, pz1, pz2, pnumx, pnumy, pnumz, pn )
  1196. FILE *infile;
  1197. real *px1, *px2, *py1, *py2, *pz1, *pz2;
  1198. int *pnumx, *pnumy, *pnumz, *pn;
  1199. #endif
  1200. {
  1201.     char tempStr[40];
  1202.     /* ----------------------------------
  1203.         Read the header information
  1204.         ---------------------------------- */
  1205.     ReadLn(infile);        /*first line is a title - ignore */
  1206.     if (fgets(tempStr, 40, infile)!= NULL)
  1207.     {
  1208.         /* Second line should have grid bounds information */
  1209.         if (6 != sscanf(tempStr, "%g %g %g %g %g %g", \
  1210.         px1, px2, py1, py2, pz1, pz2))
  1211.         {
  1212.             FileFormatError("ReadGSNFile", "poor format for grid bounds");
  1213.             return (-1);
  1214.         }
  1215.     }
  1216.     else
  1217.     {    
  1218.         FileFormatError("ReadGSNFile", "grid bounds expected");
  1219.         return (-1);
  1220.     }
  1221.     if (fgets(tempStr, 40, infile)!= NULL)
  1222.     {
  1223.         /* Third line should have grid spacing information */
  1224.         if (3 != sscanf(tempStr, "%d %d %d", pnumx, pnumy, pnumz))
  1225.         {
  1226.             FileFormatError("ReadGSNFile", "poor format for grid spacing");
  1227.             return (-1);
  1228.         }
  1229.     }
  1230.     else
  1231.     {    
  1232.         FileFormatError("ReadGSNFile", "grid spacing expected");
  1233.         return (-1);
  1234.     }
  1235.     if (fgets(tempStr, 40, infile)!= NULL)
  1236.     {
  1237.         /* Fourth line should have number of terms in the summation series */
  1238.         if (1 != sscanf(tempStr, "%d", pn))
  1239.         {
  1240.             FileFormatError("ReadGSNFile", "poor format for number of terms in series");
  1241.             return (-1);
  1242.         }
  1243.     }
  1244.     else
  1245.     {    
  1246.         FileFormatError("ReadGSNFile", "number of terms in series expected");
  1247.         return (-1);
  1248.     }
  1249. }
  1250.  
  1251. /* ----------------------------------------
  1252.     ReadGSN1Parameters
  1253.     ---------------------------------------- */    
  1254. #ifdef PROTO
  1255. static int ReadGSN1Parameters(FILE *infile, gaussianParameters *p, int n)
  1256. #else
  1257. static int ReadGSN1Parameters(infile, p, n)
  1258. FILE *infile;
  1259. gaussianParameters *p;
  1260. int n;
  1261. #endif
  1262. {
  1263.     int i;
  1264.     char tempStr[80];
  1265.  
  1266.     for (i = 0; i < n; i++ )
  1267.     {
  1268.         if (fgets (tempStr, 80, infile) == NULL )
  1269.         {
  1270.             FileFormatError("ReadGSNFile", "Expected more parameters for gaussian \
  1271. function");
  1272.             return(-1);
  1273.         }
  1274.  
  1275.  
  1276.         else  if ( 8 == sscanf(tempStr, "%g %g %d %d %d %g %g %g",  &(p+i)->c, &(p+i)->zeta, \
  1277.         &(p+i)->nx, &(p+i)->ny, &(p+i)->nz, &(p+i)->Rx,&(p+i)->Ry, &(p+i)->Rz) )
  1278.         {
  1279.             (p+i)->d = 1.0;    /* 8 column format OK */
  1280.         }
  1281.         else
  1282.         {
  1283.             FileFormatError("ReadGSNFile", "Poor format for gaussian parameters");
  1284.             return(-1);
  1285.         }
  1286.     }
  1287. }
  1288.  
  1289. /* ---------------------------------------------
  1290.     ReadGSN1File
  1291.     GSN1 file reader
  1292.  
  1293.     File format:
  1294.        **************************************************
  1295.     TITLE
  1296.     xmin xmax ymin ymax zmin zmax
  1297.     nGridx nGridy nGridz
  1298.     n
  1299.     c1 zeta1 nx1 ny1 nz1 Rx1 Ry1 Rz1
  1300.     c2 zeta2 nx2 ny2 nz2 Rx2 Ry2 Rz2
  1301.     .
  1302.       .
  1303.     .
  1304.     cn zetan nxn nyn nzn Rxn Ryn Rzn
  1305.     *************************************************
  1306.   ----------------------------------------------- */
  1307. #ifdef PROTO
  1308. static ObjPtr ReadGSN1File(ObjPtr reader, char *name)
  1309. #else
  1310. static ObjPtr ReadGSN1File(reader,name)        /* 8 - column format */
  1311. ObjPtr  reader;
  1312. char *name;
  1313. #endif
  1314. /*Reads a GSN1 file from name*/
  1315. {
  1316.     FILE *infile;
  1317.         ObjPtr retVal;                
  1318.     ObjPtr dataForm;
  1319.     ObjPtr gridDataset;
  1320.       long    dims[3];
  1321.       real    bounds[6];    /* Used to define parameters of datasets */
  1322.       long    indices[3];    /* and dataforms */
  1323.     char formName[80], dataName[80];
  1324.     char *tc;
  1325.     real curData;
  1326.     char tempStr[40];
  1327.     real min, max, minMax[2];
  1328.     ObjPtr minMaxArray;    
  1329.     real xmin, xmax, ymin, ymax, zmin, zmax;
  1330.     int ngridx, ngridy, ngridz, n;
  1331.     gaussianParameters *parameters;
  1332.     real x, y, z;
  1333.     real xInc, yInc, zInc;
  1334.     int i;
  1335.  
  1336.  
  1337.       /*---------------------
  1338.          Derive a dataset name from the file name 
  1339.        --------------------- */
  1340.  
  1341.       strcpy(formName, name);
  1342.       tc = formName;
  1343.       while (*tc && (*tc != '@') && (*tc != '.'))
  1344.         tc++;
  1345.       *tc = '\0';
  1346.       strcpy(dataName, formName);
  1347.  
  1348.       strcat(formName, ".form");
  1349.       strcat(dataName, ".data");
  1350.   
  1351.  
  1352.   /* Open the file and verify that it contains appropriate data */
  1353.  
  1354.       if ((infile = fopen(name, "r")) == NULL) 
  1355.     {
  1356.             FileFormatError("ReadGSNFile", "File not accessible");
  1357.             return NULLOBJ;
  1358.       }
  1359.  
  1360.  
  1361.     if(ParseGSN12Header(infile, &xmin, &xmax, &ymin, &ymax, &zmin, &zmax, \
  1362.     &ngridx, &ngridy, &ngridz, &n) == -1)
  1363.     {
  1364.         FileFormatError("ReadGSNFile", "Unable to parse grid parameters");
  1365.         return NULLOBJ;
  1366.     }
  1367.     /* -----------------------------------
  1368.         Allocate memory for array of function parameters and read them in
  1369.         ----------------------------------- */
  1370.     parameters = (gaussianParameters *) Alloc ( n * sizeof (gaussianParameters) );
  1371.     if (ReadGSN1Parameters(infile, parameters, n) == -1)
  1372.     {
  1373.         FileFormatError("ReadGSNFile", "Unable to read parameters");
  1374.         return NULLOBJ;
  1375.     }
  1376.  
  1377.  
  1378.     /* -----------------------------------
  1379.         Calculate the dimensions of the grid.
  1380.         There is always one more grid point than number
  1381.         of grid spaces in each direction.
  1382.        ------------------------------------ */
  1383.     dims[0] = ngridx + 1;
  1384.     dims[1] = ngridy + 1;
  1385.     dims[2] = ngridz + 1;
  1386.  
  1387.     /* -----------------------------------
  1388.         Create the scalar dataset
  1389.         ----------------------------------- */
  1390.     retVal = NewStructuredDataset(dataName, 3, dims, 0 );
  1391.  
  1392.     /* ------------------------------------
  1393.         Calculate grid bounds 
  1394.         ------------------------------------ */
  1395.     bounds[0] =  xmin;
  1396.     bounds[1] = xmax;
  1397.     bounds[2] =  ymin;
  1398.     bounds[3] = ymax;
  1399.     bounds[4] =  zmin;
  1400.     bounds[5] = zmax;
  1401.  
  1402.     /* ------------------------------------
  1403.         Create a regular data form
  1404.         ------------------------------------  */
  1405.     dataForm = NewRegularDataForm(formName, 3, dims, bounds);
  1406.     SetDatasetForm(retVal, dataForm);
  1407.     /* ------------------------------------
  1408.         Fill up the dataset (with x varying fastest)
  1409.         generating data values as you go
  1410.         Keep track of minimum and maximum values as
  1411.         data values are read in.
  1412.         ------------------------------------ */
  1413.     SetCurField(FIELD1, retVal);  /* Link the field to the data form */
  1414.     min = PLUSINF;                /* Initialize min and max to */
  1415.     max = MINUSINF;           /* extreme values */    
  1416.  
  1417.     x = xmin;
  1418.     y = ymin;
  1419.     z = zmin;
  1420.     xInc = (xmax - xmin ) / ( (real) (ngridx));
  1421.     yInc = (ymax - ymin) / ( (real) (ngridy));
  1422.     zInc = (zmax - zmin) / ((real) (ngridz));
  1423.  
  1424.     for (indices[2] = 0; indices[2] < dims[2]; (indices[2])++)
  1425.     {
  1426.         y = ymin;
  1427.         for (indices[1] = 0; indices[1] < dims[1]; (indices[1])++)
  1428.         {
  1429.             x = xmin;
  1430.             for (indices[0] = 0; indices[0] < dims[0]; (indices[0])++)
  1431.             {
  1432.                 curData = 0.0;
  1433.                 for (i=0; i < n; i++ )
  1434.                     curData +=  parameters[i].d * parameters[i].c * \
  1435.                     gaussianFunction(x, y, z, parameters[i].Rx, parameters[i].Ry,\
  1436.                      parameters[i].Rz, parameters[i].zeta, parameters[i].nx, \
  1437.                      parameters[i].ny, parameters[i].nz, 1.0);
  1438.  
  1439.             /*    printf("%f ", curData);*/
  1440.  
  1441.                 PutFieldComponent(FIELD1, 0, indices,curData);
  1442.                 if (curData < min)
  1443.                     min = curData;
  1444.                 if (curData > max)
  1445.                     max = curData;
  1446.                 x += xInc;
  1447.             }
  1448.     /*         printf("\n"); */
  1449.             y += yInc;
  1450.         }
  1451.         z += zInc;
  1452.     }
  1453.  
  1454.     /* -------------------------
  1455.         Record the minimum and maximum data values
  1456.         for colour palette and isosurface scales 
  1457.         ------------------------- */
  1458.  
  1459.     minMax[0] = min;
  1460.     minMax[1] = max; 
  1461.     minMaxArray = NewRealArray(1, 2L);
  1462.     CArray2Array(minMaxArray, minMax);
  1463.     SetVar(retVal, MINMAX, minMaxArray);
  1464.  
  1465.     RegisterDataset(retVal);
  1466.     fclose(infile);
  1467.     return(NULLOBJ);
  1468. }
  1469. /* -----------------------------------------
  1470.     ReadGSN2Parameters
  1471.     ----------------------------------------- */
  1472. #ifdef PROTO
  1473. static int ReadGSN2Parameters(FILE *infile, gaussianParameters *p, int n)
  1474. #else
  1475. static int ReadGSN2Parameters(infile, p, n)    /* 9 column */
  1476. FILE *infile;
  1477. gaussianParameters *p;
  1478. int n;
  1479. #endif
  1480. {
  1481.     int i;
  1482.     char tempStr[80];
  1483.  
  1484.     for (i = 0; i < n; i++ )
  1485.     {
  1486.         if (fgets (tempStr, 80, infile) == NULL )
  1487.         {
  1488.             FileFormatError("ReadGSNFile", "Expected more parameters for gaussian \
  1489. function");
  1490.             return(-1);
  1491.         }
  1492.  
  1493.         if (9 == sscanf(tempStr, "%g %g %g %d %d %d %g %g %g", &(p+i)->d, &(p+i)->c, &(p+i)->zeta, \
  1494.         &(p+i)->nx, &(p+i)->ny, &(p+i)->nz, &(p+i)->Rx,&(p+i)->Ry, &(p+i)->Rz) )
  1495.         { 
  1496.             /* 9 column format OK */
  1497.         }
  1498.         else
  1499.         {
  1500.             FileFormatError("ReadGSNFile", "Poor format for gaussian parameters");
  1501.             return(-1);
  1502.         }
  1503.     }
  1504. }
  1505.  
  1506. /* ---------------------------------------------
  1507.     ReadGSN2File
  1508.  
  1509.     File format:
  1510.        **************************************************
  1511.     TITLE
  1512.     xmin xmax ymin ymax zmin zmax
  1513.     nGridx nGridy nGridz
  1514.     n
  1515.     d1 c1 zeta1 nx1 ny1 nz1 Rx1 Ry1 Rz1
  1516.     d2 c2 zeta2 nx2 ny2 nz2 Rx2 Ry2 Rz2
  1517.     .
  1518.       .
  1519.     .
  1520.     dn cn zetan nxn nyn nzn Rxn Ryn Rzn
  1521.     *************************************************
  1522.   ------------------------------------------------ */
  1523. #ifdef PROTO
  1524. static ObjPtr ReadGSN2File(ObjPtr reader, char *name)
  1525. #else
  1526. static ObjPtr ReadGSN2File(reader,name)        /* 9 - column format */
  1527. ObjPtr  reader;
  1528. char *name;
  1529. #endif
  1530. /*Reads a GSN file from name*/
  1531. {
  1532.     FILE *infile;
  1533.         ObjPtr retVal;                
  1534.     ObjPtr dataForm;
  1535.     ObjPtr gridDataset;
  1536.       long    dims[3];
  1537.       real    bounds[6];    /* Used to define parameters of datasets */
  1538.       long    indices[3];    /* and dataforms */
  1539.     char formName[80], dataName[80];
  1540.     char *tc;
  1541.     real curData;
  1542.     char tempStr[40];
  1543.     real min, max, minMax[2];
  1544.     ObjPtr minMaxArray;    
  1545.     real xmin, xmax, ymin, ymax, zmin, zmax;
  1546.     int ngridx, ngridy, ngridz, n;
  1547.     gaussianParameters *parameters;
  1548.     real x, y, z;
  1549.     real xInc, yInc, zInc;
  1550.     int i;
  1551.  
  1552.  
  1553.       /*---------------------
  1554.          Derive a dataset name from the file name 
  1555.        --------------------- */
  1556.  
  1557.       strcpy(formName, name);
  1558.       tc = formName;
  1559.       while (*tc && (*tc != '@') && (*tc != '.'))
  1560.         tc++;
  1561.       *tc = '\0';
  1562.       strcpy(dataName, formName);
  1563.  
  1564.       strcat(formName, ".form");
  1565.       strcat(dataName, ".data");
  1566.   
  1567.  
  1568.   /* Open the file and verify that it contains appropriate data */
  1569.  
  1570.       if ((infile = fopen(name, "r")) == NULL) 
  1571.     {
  1572.             FileFormatError("ReadGSNFile", "File not accessible");
  1573.             return NULLOBJ;
  1574.       }
  1575.  
  1576.  
  1577.     if(ParseGSN12Header(infile, &xmin, &xmax, &ymin, &ymax, &zmin, &zmax, \
  1578.     &ngridx, &ngridy, &ngridz, &n) == -1)
  1579.     {
  1580.         FileFormatError("ReadGSNFile", "Unable to parse grid parameters");
  1581.         return NULLOBJ;
  1582.     }
  1583.     /* -----------------------------------
  1584.         Allocate memory for array of function parameters and read them in
  1585.         ----------------------------------- */
  1586.     parameters = (gaussianParameters *) Alloc ( n * sizeof (gaussianParameters) );
  1587.     if (ReadGSN2Parameters(infile, parameters, n) == -1)
  1588.     {
  1589.         FileFormatError("ReadGSNFile", "Unable to read parameters");
  1590.         return NULLOBJ;
  1591.     }
  1592.  
  1593.  
  1594.     /* -----------------------------------
  1595.         Calculate the dimensions of the grid.
  1596.         There is always one more grid point than number
  1597.         of grid spaces in each direction.
  1598.        ------------------------------------ */
  1599.     dims[0] = ngridx + 1;
  1600.     dims[1] = ngridy + 1;
  1601.     dims[2] = ngridz + 1;
  1602.  
  1603.     /* -----------------------------------
  1604.         Create the scalar dataset
  1605.         ----------------------------------- */
  1606.     retVal = NewStructuredDataset(dataName, 3, dims, 0 );
  1607.  
  1608.     /* ------------------------------------
  1609.         Calculate grid bounds 
  1610.         ------------------------------------ */
  1611.     bounds[0] =  xmin;
  1612.     bounds[1] = xmax;
  1613.     bounds[2] =  ymin;
  1614.     bounds[3] = ymax;
  1615.     bounds[4] =  zmin;
  1616.     bounds[5] = zmax;
  1617.  
  1618.     /* ------------------------------------
  1619.         Create a regular data form
  1620.         ------------------------------------  */
  1621.     dataForm = NewRegularDataForm(formName, 3, dims, bounds);
  1622.     SetDatasetForm(retVal, dataForm);
  1623.     /* ------------------------------------
  1624.         Fill up the dataset (with x varying fastest)
  1625.         generating data values as you go
  1626.         Keep track of minimum and maximum values as
  1627.         data values are read in.
  1628.         ------------------------------------ */
  1629.     SetCurField(FIELD1, retVal);  /* Link the field to the data form */
  1630.     min = PLUSINF;                /* Initialize min and max to */
  1631.     max = MINUSINF;           /* extreme values */    
  1632.  
  1633.     x = xmin;
  1634.     y = ymin;
  1635.     z = zmin;
  1636.     xInc = (xmax - xmin ) / ( (real) (ngridx));
  1637.     yInc = (ymax - ymin) / ( (real) (ngridy));
  1638.     zInc = (zmax - zmin) / ((real) (ngridz));
  1639.  
  1640.     for (indices[2] = 0; indices[2] < dims[2]; (indices[2])++)
  1641.     {
  1642.         y = ymin;
  1643.         for (indices[1] = 0; indices[1] < dims[1]; (indices[1])++)
  1644.         {
  1645.             x = xmin;
  1646.             for (indices[0] = 0; indices[0] < dims[0]; (indices[0])++)
  1647.             {
  1648.                 curData = 0.0;
  1649.                 for (i=0; i < n; i++ )
  1650.                     curData +=  parameters[i].d * parameters[i].c * \
  1651.                     gaussianFunction(x, y, z, parameters[i].Rx, parameters[i].Ry,\
  1652.                      parameters[i].Rz, parameters[i].zeta, parameters[i].nx, \
  1653.                      parameters[i].ny, parameters[i].nz, 1.0);
  1654.  
  1655.  
  1656.  
  1657.                 PutFieldComponent(FIELD1, 0, indices,curData);
  1658.                 if (curData < min)
  1659.                     min = curData;
  1660.                 if (curData > max)
  1661.                     max = curData;
  1662.                 x += xInc;
  1663.             }
  1664.  
  1665.             y += yInc;
  1666.         }
  1667.         z += zInc;
  1668.     }
  1669.  
  1670.     /* -------------------------
  1671.         Record the minimum and maximum data values
  1672.         for colour palette and isosurface scales 
  1673.         ------------------------- */
  1674.  
  1675.     minMax[0] = min;
  1676.     minMax[1] = max; 
  1677.     minMaxArray = NewRealArray(1, 2L);
  1678.     CArray2Array(minMaxArray, minMax);
  1679.     SetVar(retVal, MINMAX, minMaxArray);
  1680.  
  1681.     RegisterDataset(retVal);
  1682.     fclose(infile);
  1683.     return(NULLOBJ);
  1684. }
  1685. /* ------------------------
  1686.     FtoC()
  1687.     Converts a string rep. floating point number from FORTRAN to C
  1688.     by replacing each 'D' with an E
  1689.  --------------------------*/
  1690. #ifdef PROTO
  1691. static real FtoC(char *s1)
  1692. #else
  1693. static real FtoC(s1)
  1694. char *s1;
  1695. #endif
  1696. {
  1697.     char *tmpstr;
  1698.     real val;
  1699.  
  1700.     tmpstr = s1;
  1701.     while((*tmpstr != '\0') &&  (*tmpstr != 'D'))
  1702.         tmpstr++;
  1703.     if (*tmpstr == 'D')
  1704.         *tmpstr = 'E';
  1705.     sscanf(s1, "%f", &val);
  1706.  
  1707.     return(val);
  1708. }
  1709.  
  1710.     
  1711. /* -----------------------------------------------
  1712.     ParseGSNFile
  1713.  
  1714.     Extracts information required to generate gaussian parameters 
  1715.     from G90 output.
  1716.     
  1717.     There are six states of the parser:
  1718.     state 0:    Search for the keywords, "Standard orientation:\n"
  1719.             Read the table and store the coordinates and 
  1720.             atomic number of each atom of the molecule in
  1721.             atomArray, and the number of atoms in numAtoms.
  1722.  
  1723.     state 1:    Search for "Basis set in the form of general basis input:\n"
  1724.             What follows will be the exponents and SPD coeffs for 
  1725.             the basis functions (+diffuse and/or polarization functions)
  1726.             for each type of shell for each atom in order.
  1727.             Store the number of shell types per atom (S, SP, P, D) in 
  1728.             the array, atomShells.  Store the exponents, coefficients and
  1729.             shell type for each shell in the array, basisSet, for all atoms
  1730.             in order.  The size of basisSet[] will be the total number of
  1731.             shells (numShells), which is the sum of the number of
  1732.             shells per atom over all atoms.  The allocated size of basisSet[]
  1733.             is based on the premise that the maximum number of different
  1734.             shell types for a 3rd row atom is 6; the allocated size is then
  1735.             6 * numAtoms.
  1736.  
  1737.     state 2:    Search for "There are %d symmetry adapted basis functions..."
  1738.             There will be a variable number of lines of this form.  
  1739.             Add up the number of symmetries to determine the number 
  1740.             of molecular orbitals.
  1741.             Store this value in numSym.
  1742.  
  1743.     state 3:    Search for "%d basis functions %d primitive gaussians"
  1744.             These values are the number of orbital coefficients given
  1745.             per orbital, and the total number of terms in the summation.
  1746.  
  1747.     state 4:    Search for "Molecular Orbital Coefficients"
  1748.             which is essentially a numBasis x numSym array of coefficients,
  1749.             grouped in columns of 5 with a 2-3 line header for each set.
  1750.             Read to the orbital number desired, and then store the
  1751.             column of coefficients in orbCoeffs[].
  1752.  
  1753.     state 5:     done
  1754.  ------------------------------------------------- */             
  1755. #ifdef PROTO
  1756. static int ParseGSNFile(FILE *inFile, nuclearCentre **atomArray, int *numAtoms, \
  1757. int *numSym, int *numBasis, int *numPrimitives, basis **basisSet, int *numShells,\
  1758. int **atomShells,real **orbCoeff, int numOrb)
  1759. #else
  1760. static int ParseGSNFile(inFile, atomArray, numAtoms, numSym, numBasis, numPrimitives,\
  1761. basisSet, numShells, sPerA, orbCoeff, numOrb)
  1762. FILE *inFile;
  1763. nuclearCentre **atomArray;
  1764. int *numAtoms, *numSym, *numBasis, *numPrimitives;
  1765. basis **basisSet;
  1766. int *numShells;
  1767. int **atomShells;    /* shells per atom */
  1768. real **orbCoeff;
  1769. int numOrb;
  1770. #endif
  1771. {
  1772.     char tmpStr[80], s1[80], s2[80];
  1773.     char field[3][20];
  1774.     int i, num,done, count, count2,m, done2;
  1775.     real x, y, z;
  1776.     nuclearCentre *a;    /* pointer to array of atoms */
  1777.     int atomsAllocated = 10;    
  1778.     basis *b;        /*pointer to 1D array of basis set info */
  1779.     char *tmp;    
  1780.     int *shells;    /* pointer to number of shells per atom */
  1781.     int numSets, remainder,  j, k, lower, upper;
  1782.     real c[5];        /* temporarily holds molecular orbital coeffs */
  1783.     int state = 0;
  1784.     real *array;    /* pointer to 1D array of MOC column */
  1785.     int setDesired, columnDesired;
  1786.  
  1787.     *numSym = 0, count=0;
  1788.     SkipBlanks(inFile);
  1789.     count2=0;
  1790.     while    (fgets(tempStr, 80, inFile)!= NULL)
  1791.      {    
  1792.         
  1793.         switch (state) 
  1794.         {
  1795.             case 0:                
  1796.                 if (strcmp(tempStr, "Standard orientation:\n") == 0)
  1797.                 {        /* parse the atomic numbers and coordinates */
  1798.                         /* ------------------
  1799.                                Allocate atom array 
  1800.                            ------------------ */
  1801.                     a = (nuclearCentre *)Alloc(atomsAllocated * sizeof(nuclearCentre) );
  1802.                     if (a == NULL)
  1803.                     {
  1804.                         FileFormatError("ReadGSNFile", "Unable to allocate memory for\
  1805.  atoms.");
  1806.                         return(-1);
  1807.                     }
  1808.                     
  1809.                     ReadLn(inFile);
  1810.                     ReadLn(inFile);
  1811.                     ReadLn(inFile);
  1812.                     ReadLn(inFile);
  1813.                     done = FALSE;
  1814.                     *numAtoms = 0;
  1815.                     while((fgets(tempStr, 80, inFile) != NULL ) && !done)
  1816.                     {
  1817.                         if (sscanf(tempStr, "%d %d %g %g %g", &i, &num, &x, &y, &z)\
  1818.                         == 5)
  1819.                         {
  1820.                             /* Check if there's room in the array */
  1821.                             if ( *numAtoms  == atomsAllocated)
  1822.                             {
  1823.                                 atomsAllocated += 10;
  1824.                                 a = (nuclearCentre *)Realloc(a, atomsAllocated\
  1825.                              * sizeof(nuclearCentre) );
  1826.                                 if ( a== NULL)
  1827.                                 {
  1828.                                     FileFormatError("ReadGSNFile", \
  1829. "Could not allocate memory for atoms.");
  1830.                                     return(-1);
  1831.                                 }
  1832.                             } 
  1833.             
  1834.                             (a+i-1)->atomicNumber = num;
  1835.                             (a+i-1)->x = x;
  1836.                             (a+i-1)->y = y;    
  1837.                             (a+i-1)->z = z;
  1838.                             (*numAtoms)++;
  1839.                         }
  1840.                         else
  1841.                             done = TRUE;
  1842.                     }
  1843.                     *atomArray = a;
  1844.                     state++;
  1845.                 }
  1846.                     
  1847.                 break;
  1848.             case 1:    /* read in basis set information */
  1849.                 if (strcmp(tempStr, "Basis set in the form of general basis input:\n") == 0)
  1850.                 {
  1851.                     /* ----------------
  1852.                         Allocate memory for basis set array
  1853.                         Assume max # shells for a third row atom is 6
  1854.                         Add a safety factor of 1
  1855.                         Assume each atom is a third row (although we
  1856.                         can figure this out)
  1857.                         Also allocate for array recording shells per atom
  1858.                      ----------------- */
  1859.                     b = (basis *) Alloc ( (6+1) * (*numAtoms) * sizeof(basis) );
  1860.                     if (b == NULL)
  1861.                     {
  1862.                         FileFormatError("ReadGSNFile", \
  1863.                         "Unable to allocate memory for basis set");
  1864.                         return(-1);        /*should print error here */
  1865.                     }
  1866.                     shells = (int *) Alloc( (*numAtoms) * sizeof(int));
  1867.  
  1868.                     count = 0;    /* atom # */
  1869.                     i = 0;        /* shell # */
  1870.                     while ( count < (*numAtoms))    /*read basis set for each atom */
  1871.                     {
  1872.                         m = 0;    /* # shells for current atom */
  1873.                         ReadLn(inFile);    /* first line is centre (atom) number; discard */
  1874.                         SkipBlanks(inFile);
  1875.                         if (fgets(tempStr, 80, inFile) == NULL )
  1876.                         {
  1877.                             FileFormatError("ReadGSNFile",\
  1878.                             "Expected more data");
  1879.                             return(-1);
  1880.                         }
  1881.                         while ( strcmp(tempStr, "****\n") )
  1882.                         {    /*loop for each shell of current atom */
  1883.                             sscanf(tempStr, "%s %d %f", &(b+i)->shellType, &(b+i)->numPrim,\
  1884.                             &(b+i)->scaleFactor);
  1885.                             /* -----------------
  1886.                                 Allocate memory for exponent and coeff arrays 
  1887.                                 ----------------- */
  1888.                             b[i].exponents = (real *) Alloc ( b[i].numPrim * sizeof(real));
  1889.                             b[i].s = NULL;    /* occasionally not auto-initialized to NULL */
  1890.                             b[i].p = NULL;
  1891.                             b[i].d = NULL;
  1892.                             tmp = b[i].shellType;
  1893.                             while( *tmp != '\0')
  1894.                             {
  1895.                                 switch(*tmp)
  1896.                                 {
  1897.                                     case 'S':
  1898.                                         b[i].s = (real *) Alloc ( b[i].numPrim * sizeof(real));
  1899.                                         break;
  1900.                                     case 'P':
  1901.                                         b[i].p = (real *) Alloc ( b[i].numPrim * sizeof(real));
  1902.                                         break;
  1903.                                     case 'D':
  1904.                                         b[i].d = (real *) Alloc ( b[i].numPrim * sizeof(real));
  1905.                                         break;
  1906.                                     default:
  1907.                                         FileFormatError("ReadGSNFile", \
  1908.                                         "Unknown shell type");
  1909.                                         return(-1);
  1910.                                         break;
  1911.                                 }
  1912.                                 tmp++;
  1913.                             }
  1914.  
  1915.                             for (j = 0; j < b[i].numPrim; j++ )
  1916.                             {
  1917.                                 if (fgets(tempStr, 80, inFile) == NULL )
  1918.                                 {
  1919.                                     FileFormatError("ReadGSNFile",\
  1920.                                     "Expected more data");
  1921.                                     return(-1);
  1922.                                 }
  1923.                                 sscanf(tempStr, "%s %s %s", field[0], field[1], field[2]);
  1924.                                 b[i].exponents[j] = FtoC(field[0]);
  1925.                             
  1926.                                 k = 1; /*field array index*/
  1927.                                 tmp = b[i].shellType;
  1928.                                 while( *tmp != '\0')
  1929.                                 {
  1930.                                     switch(*tmp)
  1931.                                     {
  1932.                                         case 'S':
  1933.                                             b[i].s[j] = FtoC(field[k]); k++;
  1934.                                             break;
  1935.                                         case 'P':
  1936.                                             b[i].p[j] = FtoC(field[k]); k++;
  1937.                                             break;
  1938.                                         case 'D':
  1939.                                             b[i].d[j] = FtoC(field[k]); k++;
  1940.                                             break;
  1941.                                         default:
  1942.                                             FileFormatError("ReadGSNFile", \
  1943.                                             "Unknown shell type");
  1944.                                             return(-1);
  1945.                                             break;    
  1946.                                     }
  1947.                                     tmp++;
  1948.                                 } /* end while (for cycling through shell type, fields)*/
  1949.                             } /* end for (for reading exp, coeff for each primitive)*/
  1950.                             i++;    /* increment shell # */
  1951.                             m++;
  1952.                             SkipBlanks(inFile);
  1953.                             if (fgets(tempStr, 80, inFile) == NULL )
  1954.                             {
  1955.                                 FileFormatError("ReadGSNFile",\
  1956.                                 "Expected more data");
  1957.                                 return(-1);
  1958.                             }
  1959.                         } /* end while (for reading each shell for current atom) */
  1960.                         shells[count]= m;
  1961.                         count++;
  1962.                     } /* end while (for reading basis sets for each atom ) */
  1963.                     *basisSet = b;
  1964.                     *numShells = i;
  1965.                     *atomShells = shells;
  1966.                     
  1967.                     state++;
  1968.                 } /* end if (for checking for basis set keywords ) */
  1969.                 break;
  1970.             case 2:
  1971.                 if (sscanf(tempStr, \
  1972.                 "There are %d symmetry adapted basis functions of %s symmetry.",\
  1973.                  &i, s1) == 2 )
  1974.                 {    /* increment the count of number of orbital symmetries 
  1975.                         while the "There are..." pattern is matched */
  1976.                     *numSym = i;
  1977.  
  1978.                     SkipBlanks(inFile);
  1979.                     if(fgets(tempStr, 80, inFile)== NULL)
  1980.                     {
  1981.                         FileFormatError("ReadGSNFile", "Expected more data");
  1982.                         return(-1);
  1983.                     }
  1984.                     while(sscanf(tempStr, \
  1985.                     "There are %d symmetry adapted basis functions of %s symmetry.",\
  1986.                      &i, s1) == 2 )
  1987.                     {
  1988.                         (*numSym) += i;    
  1989.                         SkipBlanks(inFile);
  1990.                         if(fgets(tempStr, 80, inFile)== NULL)
  1991.                         {
  1992.                             FileFormatError("ReadGSNFile", "Expected more data");
  1993.                             return(-1);
  1994.                         }
  1995.                     }
  1996.  
  1997.                     state++;
  1998.                 }
  1999.                 
  2000.                 break;
  2001.             case 3:
  2002.                 if (sscanf(tempStr, \
  2003.                 "%d basis functions %d primitive gaussians", numBasis, \
  2004.                 numPrimitives) == 2)
  2005.                     state++;
  2006.                 break;
  2007.             case 4:
  2008.                 if (strcmp (tempStr, "Molecular Orbital Coefficients\n") == 0)
  2009.                 {
  2010.  
  2011.                     numSets = (*numSym) / 5;        /* columns are printed in sets of 5 */
  2012.                     remainder = (*numSym) % 5;    /* last set will have "remainder"# of columns */
  2013.  
  2014.                     /* ------------------
  2015.                         Check if orbital desired falls within range
  2016.                              ------------------ */
  2017.  
  2018.                     if ( (numOrb < 0) || (numOrb > (*numSym-1)) )
  2019.                     {
  2020.                         FileFormatError("ReadGSNFile", "Invalid orbital selected");
  2021.                         return(-1);
  2022.                     }
  2023.  
  2024.                     array= (real *) Alloc ((*numBasis) * sizeof(real));
  2025.                     if (array == NULL)
  2026.                         return(-1);
  2027.  
  2028.                     setDesired = numOrb  / 5;            /* 0..*numBasis-1*/
  2029.                     columnDesired = numOrb % 5;         /* 0..4 */
  2030.  
  2031.                     /* Read past all sets until set desired */        
  2032.                     for ( k = 0; k < setDesired; k++)
  2033.                     {
  2034.                         /* first 2 or 3 lines of each set are header rows */
  2035.                         done2 = FALSE;
  2036.                         while (!done2 && fgets(tempStr, 80, inFile)!= NULL)
  2037.                         {
  2038.                             sscanf(tempStr, "%s", s1);
  2039.                             if (!strcmp(s1, "EIGENVALUES") )
  2040.                                 done2 = TRUE;
  2041.                         }
  2042.                         for ( m = 0; m < (*numBasis); m++)
  2043.                             ReadLn(inFile);
  2044.                     }
  2045.                     /* now should be at correct set */
  2046.                     /* first 2 or 3 lines of each set are header rows */
  2047.                     done2 = FALSE;
  2048.                     while (!done2 && fgets(tempStr, 80, inFile)!= NULL)
  2049.                     {
  2050.                         sscanf(tempStr, "%s", s1);
  2051.                         if (!strcmp(s1, "EIGENVALUES") )
  2052.                             done2 = TRUE;
  2053.                     }
  2054.                     for ( m = 0; m < (*numBasis); m++)
  2055.                     {
  2056.                         if (fgets(tempStr, 80, inFile)!= NULL)
  2057.                         {
  2058.                             strcpy(s1, tempStr+20);
  2059.                             sscanf(s1, "%g %g %g %g %g", \
  2060.                             c, c+1, c+2, c+3, c+4 );
  2061.                             array[m] = c[columnDesired];
  2062.                         }
  2063.                     }
  2064.                     *orbCoeff = array;
  2065.                     state++;
  2066.  
  2067.                 }
  2068.                 break;
  2069.  
  2070.         }
  2071.         SkipBlanks(inFile);
  2072.     }
  2073.     if (state > 4)
  2074.         return(0);
  2075.     else
  2076.         return(-1);
  2077. }
  2078.  
  2079. /* --------------------------------------------------
  2080.     makeGSNParameters
  2081.  
  2082.     Generate parameters for all primitives in summation.  Pair up the
  2083.     orbital exponents and SPD coefficients with the molecular orbital coefficients,
  2084.     generate nx, ny, nz depending on shell type, and Rx, Ry, Rz will be
  2085.     the coordinates of the centre of the current atom.  (Easiest way to understand
  2086.     algorithm is to generate a few primitives by hand).
  2087.   ----------------------------------------------------*/
  2088.  
  2089. #ifdef PROTO
  2090. static int makeGSNParameters( gaussianParameters **param,nuclearCentre *a, int numAtoms, \
  2091. int numSym, int numBasis, int numPrimitives, basis *basisSet, int numShells,\
  2092. int *atomShells,real *orbCoeff)
  2093. #else
  2094. static int makeGSNParameters( param,a, numAtoms, numSym, numBasis, numPrimitives,\
  2095. basisSet, numShells, sPerA, orbCoeff)
  2096. gaussianParameters **param;
  2097. nuclearCentre *a;
  2098. int numAtoms, numSym, numBasis, numPrimitives;
  2099. basis *basisSet;
  2100. int numShells;
  2101. int *atomShells;    /* shells per atom */
  2102. real *orbCoeff;
  2103. #endif
  2104. {
  2105.     gaussianParameters *p;    /* pointer to array of parameters */
  2106.     int i, j, k, count;
  2107.     int mocCount;                /* row index of molecular orbital coeffs */
  2108.     int paramCount;                /* index of parameter array */
  2109.     int atomCount, shellsPerAtomCount;
  2110.     char *tmp;
  2111.     real pxMOC, pyMOC, pzMOC;
  2112.     real dxxMOC, dyyMOC, dzzMOC, dxyMOC, dxzMOC, dyzMOC;
  2113.  
  2114.     p = (gaussianParameters *) Alloc ( numPrimitives * sizeof (gaussianParameters));
  2115.     if (p == NULL)
  2116.     {
  2117.         FileFormatError ("ReadGSNFile", "Unable to allocate memory for parameters");
  2118.         return(-1);
  2119.     }
  2120.  
  2121.     /* ------------------------
  2122.         Fill in array
  2123.       ------------------------- */
  2124.     mocCount = 0;
  2125.     paramCount = 0;
  2126.     atomCount = 0;
  2127.     shellsPerAtomCount = 0;
  2128.     for (i = 0; i < numShells; i++)
  2129.     {    /* loop through basis set in general basis set input format */
  2130.         tmp = basisSet[i].shellType;
  2131.         while( *tmp != '\0')
  2132.         {
  2133.             switch(*tmp)
  2134.             {
  2135.                 case 'S':    /* get one molecular orbital coeff (constant)*/
  2136.                     /* loop through #primitives for this S shell */
  2137.                     for (j = 0; j < basisSet[i].numPrim; j++)
  2138.                     {
  2139.                         p[paramCount].d = orbCoeff[mocCount];
  2140.                         p[paramCount].c = basisSet[i].s[j];
  2141.                         p[paramCount].zeta = basisSet[i].exponents[j] * basisSet[i].scaleFactor;
  2142.                         p[paramCount].nx = 0;
  2143.                         p[paramCount].ny = 0;
  2144.                         p[paramCount].nz = 0;
  2145.  
  2146.                         p[paramCount].Rx = a[atomCount].x;
  2147.                         p[paramCount].Ry = a[atomCount].y;
  2148.                         p[paramCount].Rz = a[atomCount].z;
  2149.  
  2150.                         paramCount++;
  2151.                     }
  2152.                     mocCount++;        
  2153.                     break;
  2154.                 case 'P': /* for each primitive, get Px, Py, Pz MOC */
  2155.                     pxMOC = orbCoeff[mocCount];
  2156.                     mocCount++;
  2157.                     pyMOC = orbCoeff[mocCount];
  2158.                     mocCount++;
  2159.                     pzMOC = orbCoeff[mocCount];
  2160.                     mocCount++;
  2161.  
  2162.                     /* record where current P set begins */
  2163.                     count = paramCount;    
  2164.                     /* loop through #primitives for  Px, Py,Pz shells */
  2165.                     for ( k = 0; k < 3; k++)
  2166.                     {
  2167.                         for (j = 0; j <  basisSet[i].numPrim; j++)
  2168.                         {    
  2169.                             p[paramCount].c = basisSet[i].p[j];
  2170.                             p[paramCount].zeta = basisSet[i].exponents[j] * basisSet[i].scaleFactor;
  2171.                             p[paramCount].nx = 0;    /*initialize to zero*/
  2172.                             p[paramCount].ny = 0;        /*change nx ny nz */
  2173.                             p[paramCount].nz = 0;        /* later */
  2174.  
  2175.                             p[paramCount].Rx = a[atomCount].x;
  2176.                             p[paramCount].Ry = a[atomCount].y;
  2177.                             p[paramCount].Rz = a[atomCount].z;
  2178.  
  2179.                             paramCount++;
  2180.                         }
  2181.                     }
  2182.                     /* loop through #primitives for  Px shell */
  2183.                     /* start at "count" */
  2184.                     for (j = 0; j < basisSet[i].numPrim; j++)
  2185.                     {
  2186.                         p[count].d = pxMOC;
  2187.                         p[count].nx = 1;
  2188.                         count++;
  2189.                     }
  2190.                     /* loop through #primitives for  Py shell */
  2191.                     for (j = 0; j < basisSet[i].numPrim; j++)
  2192.                     {
  2193.                         p[count].d = pyMOC;
  2194.                         p[count].ny = 1;
  2195.                         count++;
  2196.                     }
  2197.                     /* loop through #primitives for  Pz shell */
  2198.                     for (j = 0; j < basisSet[i].numPrim; j++)
  2199.                     {
  2200.                         p[count].d = pzMOC;
  2201.                         p[count].nz = 1;
  2202.                         count++;
  2203.                     }
  2204.  
  2205.                     break;
  2206.                 case 'D':    /* read molecular orb. coeffs for d type shells*/
  2207.                     dxxMOC = orbCoeff[mocCount];
  2208.                     mocCount++;
  2209.                     dyyMOC = orbCoeff[mocCount];
  2210.                     mocCount++;
  2211.                     dzzMOC = orbCoeff[mocCount];
  2212.                     mocCount++;
  2213.                     dxyMOC = orbCoeff[mocCount];
  2214.                     mocCount++;
  2215.                     dxzMOC = orbCoeff[mocCount];
  2216.                     mocCount++;
  2217.                     dyzMOC = orbCoeff[mocCount];
  2218.                     mocCount++;
  2219.  
  2220.                     /* record where current D set begins */
  2221.                     count = paramCount;    
  2222.                     /* loop through #primitives for all 6 D type shells */
  2223.                     for ( k = 0; k < 6; k++)
  2224.                     {
  2225.                         for (j = 0; j <  basisSet[i].numPrim; j++)
  2226.                         {    
  2227.                             p[paramCount].c = basisSet[i].d[j];
  2228.                             p[paramCount].zeta = basisSet[i].exponents[j] * basisSet[i].scaleFactor;
  2229.                             p[paramCount].nx = 0;    /*initialize to zero*/
  2230.                             p[paramCount].ny = 0;        /*change nx ny nz */
  2231.                             p[paramCount].nz = 0;        /* later */
  2232.  
  2233.                             p[paramCount].Rx = a[atomCount].x;
  2234.                             p[paramCount].Ry = a[atomCount].y;
  2235.                             p[paramCount].Rz = a[atomCount].z;
  2236.  
  2237.                             paramCount++;
  2238.                         }
  2239.                     }
  2240.                     /* loop through #primitives for  Dxx shell */
  2241.                     /* fill in d, nx,ny,nz; start at "count" */
  2242.                     for (j = 0; j < basisSet[i].numPrim; j++)
  2243.                     {
  2244.                         p[count].d = dxxMOC;
  2245.                         p[count].nx = 2;
  2246.                         count++;
  2247.                     }
  2248.                     /* loop through #primitives for  Dyy shell */
  2249.                     for (j = 0; j < basisSet[i].numPrim; j++)
  2250.                     {
  2251.                         p[count].d = dyyMOC;
  2252.                         p[count].ny = 2;
  2253.                         count++;
  2254.                     }
  2255.                     /* loop through #primitives for  Dzz shell */
  2256.                     for (j = 0; j < basisSet[i].numPrim; j++)
  2257.                     {
  2258.                         p[count].d = dzzMOC;
  2259.                         p[count].nz = 2;
  2260.                         count++;
  2261.                     }
  2262.                     /* loop through #primitives for  Dxy shell */
  2263.                     /* fill in d, nx,ny,nz; start at "count" */
  2264.                     for (j = 0; j < basisSet[i].numPrim; j++)
  2265.                     {
  2266.                         p[count].d = dxyMOC;
  2267.                         p[count].nx = 1;
  2268.                         p[count].ny = 1;
  2269.                         count++;
  2270.                     }
  2271.                     /* loop through #primitives for  Dxz shell */
  2272.                     for (j = 0; j < basisSet[i].numPrim; j++)
  2273.                     {
  2274.                         p[count].d = dxzMOC;
  2275.                         p[count].nx = 1;
  2276.                         p[count].nz = 1;
  2277.                         count++;
  2278.                     }
  2279.                     /* loop through #primitives for  Dyz shell */
  2280.                     for (j = 0; j < basisSet[i].numPrim; j++)
  2281.                     {
  2282.                         p[count].d = dyzMOC;
  2283.                         p[count].ny = 1;
  2284.                         p[count].nz = 1;
  2285.                         count++;
  2286.                     }
  2287.                     
  2288.                     break;
  2289.                 default:
  2290.                     FileFormatError("ReadGSNFile", \
  2291.                     "Unknown shell type");
  2292.                     return(-1);
  2293.                     break;
  2294.             }
  2295.             tmp++;
  2296.  
  2297.  
  2298.         } /* end while (cycling through shell types found in current shell ) */    
  2299.         shellsPerAtomCount++;
  2300.         if (shellsPerAtomCount == atomShells[atomCount])
  2301.         {
  2302.             atomCount++;
  2303.             shellsPerAtomCount = 0;
  2304.         }    
  2305.     } /* end for (looping through basis set in general basis set input format) */ 
  2306.     *param = p;
  2307.     return(0);
  2308. }
  2309.  
  2310. /* --------------------------------------------
  2311.     makeXYZDataset
  2312.  
  2313.     Generate unstructured dataset to visualize molecule in
  2314.     ball-and-or-stick form.  Taken from XYZ file reader.
  2315.     (time-independent)
  2316.   --------------------------------------------- */
  2317. #ifdef PROTO 
  2318. static void makeXYZDataset ( char *name, nuclearCentre *a, int n)
  2319. #else
  2320. static void makeXYZDataset ( name,a, n)
  2321. char *name;
  2322. nuclearCentre *a;
  2323. int n;
  2324. #endif
  2325. {  
  2326.     long numAtoms;        /* number of atoms per frame */
  2327.     long nBonds = 0;
  2328.     TwoReals *bondsBuffer = 0;
  2329.     long nBondsAllocated = 0;
  2330.     long dims[2];
  2331.  
  2332.     char s[4];
  2333.     char extension[5];
  2334.  
  2335.  
  2336.     long i, j,k, count;            /* counters*/
  2337.  
  2338.     Bool convertToAngstroms;    /* false if already in angstroms, true otherwise */    
  2339.  
  2340.     /*-------------------------------------
  2341.        variables for "balls" visualization 
  2342.         -------------------------------------*/
  2343.     /* NOTE: important to have new names for each dataset and form */
  2344.     char atomsName[80], radiiName[80],  tmpName[80],formName[80];
  2345.     char *tc;        /* terminating character */
  2346.     real bounds[6];
  2347.     char atomNameBuf[4];
  2348.     ObjPtr formVectors, dataForm, atomicNumberDataset, radiusDataset;
  2349.     real ballScaleFactor;
  2350.  
  2351.     /* -------------------------------------
  2352.        variables for bonds
  2353.        ------------------------------------- */
  2354.     real maxStickLength;
  2355.     real  sum, temp;
  2356.     real covalentRadius1, covalentRadius2, scaleFactor;
  2357.     int bondCalcType;    /* method used to calculate minimum interatomic dist for bond*/
  2358.  
  2359.     numAtoms = n;
  2360.  
  2361.         bondCalcType = CALC_COVALENTRADII;
  2362.  
  2363.         scaleFactor = FUDGEFACTOR;
  2364.  
  2365.         ballScaleFactor = 1.0;
  2366.     
  2367.     convertToAngstroms =  false;
  2368.  
  2369.  
  2370.  
  2371.  
  2372.     /* -------------------------------------
  2373.         make names for datasets 
  2374.         ------------------------------------- */
  2375.     strcpy (tmpName, name);
  2376.  
  2377.     strcpy(atomsName, tmpName);
  2378.     strcpy(radiiName, tmpName);
  2379.     strcpy(formName, tmpName);
  2380.  
  2381.     strcat(atomsName, ".atoms");
  2382.     strcat(radiiName, ".radii");
  2383.     strcat(formName, ".form");    
  2384.  
  2385.     k=1;
  2386.     strcpy(tmpName, atomsName);
  2387.     while (FindDatasetByName(tmpName))
  2388.     {
  2389.         ++k;
  2390.         sprintf(tmpName, "%s(%d)", atomsName,k);
  2391.     }
  2392.         
  2393.     if (k > 1)
  2394.     {
  2395.         sprintf(extension, "(%d)", k);
  2396.         strcat(atomsName, extension);
  2397.         strcat(radiiName, extension);
  2398.         strcat(formName, extension);
  2399.     }
  2400.  
  2401.  
  2402.  
  2403.     /* ----------------------------------------
  2404.           make bounds impossible
  2405.         ---------------------------------------- */
  2406.     bounds[0] = bounds[2] = bounds[4] = PLUSINF;
  2407.     bounds[1] = bounds[3] = bounds[5] = MINUSINF;
  2408.  
  2409.  /* Loop through atoms to revise bounds */
  2410.         for (i = 0; i < n; i++ )
  2411.         {    
  2412.             if (a[i].x < bounds[0])
  2413.                 bounds[0] = a[i].x;
  2414.             if (a[i].x > bounds[1] )
  2415.                 bounds[1] = a[i].x;
  2416.             if (a[i].y < bounds[2])
  2417.                 bounds[2] = a[i].y;
  2418.             if (a[i].y > bounds[3])
  2419.                 bounds[3] = a[i].y;
  2420.             if (a[i].z < bounds[4])
  2421.                 bounds[4] = a[i].z;
  2422.             if (a[i].z > bounds[5])
  2423.                 bounds[5] = a[i].z;
  2424.         }            
  2425.  
  2426.         nBonds = 0;
  2427.         nBondsAllocated = 200;
  2428.         bondsBuffer = (TwoReals *) Alloc(sizeof(TwoReals) * nBondsAllocated);    
  2429.  
  2430.         for ( i = 0; i < (n-1); i++ )
  2431.         {
  2432.             covalentRadius1 = atomInfo[a[i].atomicNumber-1].covalentRadius;    
  2433.             /* --------------------------------
  2434.                 If either atom is missing the covalent radius data,
  2435.                     do not draw a bond between them.
  2436.             -------------------------------- */
  2437.  
  2438.             if (covalentRadius1 != MISRAD)
  2439.             {
  2440.                 for (j = (i+1); j < numAtoms; j++ )
  2441.                 {
  2442.                     covalentRadius2 = atomInfo[a[j].atomicNumber-1].covalentRadius;
  2443.  
  2444.                     if ( covalentRadius2 != MISRAD)
  2445.                     {
  2446.                         sum = 0.0;
  2447.                         temp =  a[j].x - a[i].x;
  2448.                         sum += temp *temp;
  2449.                         temp = a[j].y - a[i].y;
  2450.                         sum += temp * temp;
  2451.                         temp = a[j].z - a[i].z;
  2452.                         sum += temp * temp;
  2453.                         if (sum <= SQUARE ( (covalentRadius1 + covalentRadius2)*scaleFactor))
  2454.                         {
  2455.                                if (nBonds >= nBondsAllocated)
  2456.                                 {
  2457.                                 nBondsAllocated += 200;
  2458.                                 bondsBuffer = (TwoReals *) Realloc(bondsBuffer,
  2459.                                 sizeof(TwoReals) * nBondsAllocated);
  2460.                                 }
  2461.                                 bondsBuffer[nBonds][0] = i;
  2462.                                 bondsBuffer[nBonds][1] = j;
  2463.                                 ++nBonds;
  2464.                         }    
  2465.                     }
  2466.                 }            
  2467.             }
  2468.         }
  2469.  
  2470.  
  2471.         /* make datasets */
  2472.         atomicNumberDataset = NewDataset( atomsName, 1, &numAtoms, 0);
  2473.  
  2474.         /* this variable signals the palette creation later */
  2475.         SetVar(atomicNumberDataset, UNITSNAME, NewString("Atomic number"));
  2476.         
  2477.         radiusDataset = NewDataset( radiiName, 1, &numAtoms, 0 );
  2478.  
  2479.         /* --------------------------------------
  2480.             make a vector dataset for the form
  2481.             -------------------------------------- */
  2482.         formVectors = NewDataset(formName, 1, &numAtoms, 3);
  2483.         dims[0] = numAtoms;
  2484.         dims[1] = nBonds;
  2485.         if (nBonds > 0)
  2486.         { /* one dimensional dataset (points and lines)  */
  2487.             dataForm = NewUnstructuredDataForm(tmpName, 1, dims, bounds,
  2488. formVectors);
  2489.         }
  2490.         else    
  2491.         {/* zero dimensional dataset (points only)*/
  2492.             dataForm = NewUnstructuredDataForm(tmpName, 0, dims, bounds,
  2493. formVectors);
  2494.         }
  2495.  
  2496.         SetDatasetForm(atomicNumberDataset, dataForm);
  2497.         SetDatasetForm(radiusDataset, dataForm);
  2498.         /* Fill it all up */
  2499.         if (nBonds)
  2500.         {
  2501.         /*Set the bonds*/
  2502.         for (k = 0; k < nBonds; ++k)
  2503.         {
  2504.             Put1Edge(dataForm, k,
  2505.             (long) bondsBuffer[k][0], (long) bondsBuffer[k][1]);
  2506.         }
  2507.         }
  2508.         SetCurField(FORMFIELD, formVectors);
  2509.         SetCurField(FIELD1, atomicNumberDataset);
  2510.         SetCurField(FIELD2, radiusDataset);
  2511.         
  2512.         for (i=0; i < numAtoms; i++ )
  2513.         {
  2514.                 PutFieldComponent(FORMFIELD, 0, &i, a[i].x);
  2515.                 PutFieldComponent(FORMFIELD, 1, &i, a[i].y);
  2516.                 PutFieldComponent(FORMFIELD, 2, &i, a[i].z);
  2517.     
  2518.             
  2519.  
  2520.                 PutFieldComponent(FIELD1, 0, &i, a[i].atomicNumber ? (real)
  2521. a[i].atomicNumber : missingData);
  2522.                 PutFieldComponent(FIELD2, 0, &i, a[i].atomicNumber ?
  2523. atomInfo[a[i].atomicNumber-1].radius : missingData );
  2524.         }
  2525.  
  2526.         /* -------------------------
  2527.             create new variables of SIZEOBJ ID
  2528.             allows sizing and scaling of atomic radii
  2529.            -------------------------- */
  2530.         SetVar(atomicNumberDataset, SIZEOBJ, radiusDataset);
  2531.         SetVar(radiusDataset, SIZEOBJ, radiusDataset);
  2532.  
  2533.         SetVar(atomicNumberDataset, POINTSTYLE, NewInt(BS_SPHERES));
  2534.         RegisterDataset ( atomicNumberDataset);
  2535.         RegisterDataset ( radiusDataset );
  2536.         RegisterDataset ( formVectors );
  2537.  
  2538.  
  2539.  
  2540.     /* ---------------------------------
  2541.         HACK!
  2542.  
  2543.         Sets default ball visualization to a sphere instead
  2544.         of an octohedron.  SetVar of POINTSTYLE of 
  2545.         dataset itself doesn't seem to work.
  2546.     ----------------------------------- */
  2547.       SetVar(visBalls, POINTSTYLE, NewInt(BS_SPHERES));
  2548.     /* ---------------------------------
  2549.         MORE HACK!
  2550.  
  2551.         Set scale factor of balls' size object
  2552.         ---------------------------------- */
  2553.     SetVar(visBalls, SIZEFACTOR, NewReal(ballScaleFactor));
  2554.  
  2555.         
  2556. }
  2557. /* ---------------------------
  2558.     Default GSN values
  2559.     --------------------------- */
  2560. #define GXMIN     -5.0    /*default grid bounds */
  2561. #define GXMAX    5.0
  2562. #define GYMIN     -5.0
  2563. #define GYMAX    5.0
  2564. #define GZMIN    -5.0
  2565. #define GZMAX    5.0
  2566. #define GGRIDX    10        /* default grid spacing */
  2567. #define GGRIDY    10
  2568. #define GGRIDZ    10
  2569. #define MOLECULARORBITAL    1    /* default orbital selected is first one */
  2570.  
  2571. /* -------------------------------
  2572.     GSN File Reader
  2573.     ------------------------------- */
  2574. #ifdef PROTO
  2575. static ObjPtr ReadGSNFile(ObjPtr reader, char *name)
  2576. #else
  2577. static ObjPtr ReadGSNFile(reader,name)
  2578. ObjPtr  reader;
  2579. char *name;
  2580. #endif
  2581. /*Reads a GSN file from name*/
  2582. {
  2583.     FILE *infile;
  2584.         ObjPtr retVal;                
  2585.     ObjPtr dataForm;
  2586.     ObjPtr gridDataset;
  2587.     ObjPtr var;
  2588.       long    dims[3];
  2589.       real    bounds[6];    /* Used to define parameters of datasets */
  2590.       long    indices[3];    /* and dataforms */
  2591.     char formName[80], dataName[80], tmpName[80];
  2592.     char *tc;
  2593.     real curData;
  2594.     char tempStr[40];
  2595.     real min, max, minMax[2];
  2596.     ObjPtr minMaxArray;    
  2597.     real xmin, xmax, ymin, ymax, zmin, zmax;
  2598.     int ngridx, ngridy, ngridz, n;
  2599.     gaussianParameters *parameters;
  2600.     real x, y, z;
  2601.     real xInc, yInc, zInc;
  2602.     int i;
  2603.     /* ----------------------
  2604.         new variables
  2605.        ---------------------- */
  2606.     nuclearCentre *atoms;        /* pointer to array of atom type and coordinate info */
  2607.     int nAtoms, nSym;        /* number of atoms and orbital symmetries */
  2608.     basis *genBasis;            /* points to array of general basis set */
  2609.     int nGen;                /* length of genBasis array (total #shells) */
  2610.     int *shellsPerAtom;        /* points to array that records #shells per atom */
  2611.     int nBasis, nPrim;            /* number of basis functions and primitive gaussians */
  2612.     real *d;                    /* nBasis x 1 array of selected molecular orbital coefficients */
  2613.     int j;
  2614.     int orbitalNumber = 1;        /* number of orbital to generate (0..nSym-1)*/
  2615.     Bool normalize;            /* true if MOC's to be normalized, false otherwise */
  2616.     real sum, N;                /* normalization coefficient */
  2617.       /*---------------------
  2618.          Derive a dataset name from the file name 
  2619.        --------------------- */
  2620.  
  2621.       strcpy(tmpName, name);
  2622.       tc = tmpName;
  2623.       while (*tc && (*tc != '@') && (*tc != '.'))
  2624.         tc++;
  2625.       *tc = '\0';
  2626.       strcpy(dataName, tmpName);
  2627.     strcpy(formName, tmpName);
  2628.  
  2629.       strcat(formName, ".form");
  2630.       strcat(dataName, ".data");
  2631.       /* -------------------------------------
  2632.       get the molecular orbital selected
  2633.        ------------------------------------- */
  2634.     var = GetIntVar("ReadGSNFile", reader, USER7);
  2635.       if (!var)
  2636.             orbitalNumber = MOLECULARORBITAL;
  2637.       else
  2638.             orbitalNumber = GetInt(var);
  2639.     /* -------------------------------------
  2640.         get state of normalization checkbox
  2641.         ------------------------------------- */
  2642.     normalize = GetPredicate(reader, USER8) ? true : false;
  2643.     /* -------------------------------------
  2644.          get grid bounds and spacing
  2645.        -------------------------------------- */
  2646.     var = GetRealVar("ReadGSNFile", reader, GSNXMIN);
  2647.     if (!var)
  2648.         xmin = GXMIN;
  2649.     else
  2650.         xmin = GetReal(var);
  2651.     var = GetRealVar("ReadGSNFile", reader, GSNXMAX);
  2652.     if (!var)
  2653.         xmax = GXMAX;
  2654.     else
  2655.         xmax = GetReal(var);
  2656.     var = GetRealVar("ReadGSNFile", reader, GSNYMIN);
  2657.     if (!var)
  2658.         ymin = GYMIN;
  2659.     else
  2660.         ymin = GetReal(var);
  2661.     var = GetRealVar("ReadGSNFile", reader, GSNYMAX);
  2662.     if (!var)
  2663.         ymax = GYMAX;
  2664.     else
  2665.         ymax = GetReal(var);
  2666.     var = GetRealVar("ReadGSNFile", reader, GSNZMIN);
  2667.     if (!var)
  2668.         zmin = GZMIN;
  2669.     else
  2670.         zmin = GetReal(var);
  2671.     var = GetRealVar("ReadGSNFile", reader, GSNZMAX);
  2672.     if (!var)
  2673.         zmax = GZMAX;
  2674.     else
  2675.         zmax = GetReal(var);
  2676.     var = GetIntVar("ReadGSNFile", reader, GSNGRIDX);
  2677.     if (!var)
  2678.         ngridx = GGRIDX;
  2679.     else
  2680.         ngridx = GetInt(var);
  2681.     var = GetIntVar("ReadGSNFile", reader, GSNGRIDY);
  2682.     if (!var)
  2683.         ngridy = GGRIDY;
  2684.     else
  2685.         ngridy = GetInt(var);
  2686.     var = GetIntVar("ReadGSNFile", reader, GSNGRIDZ);
  2687.     if (!var)
  2688.         ngridz = GGRIDZ;
  2689.     else
  2690.         ngridz = GetInt(var);
  2691. /* Check if bounds are OK */
  2692.   if ( (xmax < xmin) || (ymax < ymin) || (zmax < xmin) )
  2693.   {
  2694.     FileFormatError("ReadGSNFile", "Invalid grid bounds");
  2695.     return NULLOBJ;
  2696.    }
  2697.  
  2698.  
  2699.   /* Open the file and verify that it contains appropriate data */
  2700.  
  2701.       if ((infile = fopen(name, "r")) == NULL) 
  2702.     {
  2703.             FileFormatError("ReadGSNFile", "File not accessible");
  2704.             return NULLOBJ;
  2705.       }
  2706.  
  2707.  
  2708.     if (ParseGSNFile(infile, &atoms, &nAtoms, &nSym, &nBasis, &nPrim,\
  2709.     &genBasis, &nGen, &shellsPerAtom, &d, orbitalNumber-1) == -1 )
  2710.     {
  2711.         FileFormatError("ReadGSNFile", "Unable to parse gaussian parameters");
  2712.         return NULLOBJ;
  2713.     }
  2714.     for (i=0; i < nAtoms; i++ )
  2715.     {
  2716.         printf("%d  %d %g %g %g\n", i+1, atoms[i].atomicNumber, atoms[i].x, \
  2717.         atoms[i].y, atoms[i].z );
  2718.         
  2719.     }
  2720.     printf("Number of orbital symmetries: %d\n", nSym);
  2721.     printf("Number of basis functions: %d\n", nBasis);
  2722.     printf("Number of primitive gaussians: %d\n", nPrim);
  2723.     
  2724.     printf("Total number of shells: %d\n", nGen);
  2725.     printf("Shells per atom: \n");
  2726.     for (i = 0; i < nAtoms; i++)
  2727.         printf("Atom:%d Shells:%d\n", i+1, shellsPerAtom[i]);
  2728.  
  2729.     for (i = 0; i < nGen; i++)
  2730.     {
  2731.         printf("%d Shell type:%s Scale factor: %f Number of primitives: %d\n",\
  2732.         i, genBasis[i].shellType, genBasis[i].scaleFactor, genBasis[i].numPrim);
  2733.         for (j=0; j < genBasis[i].numPrim; j++ )
  2734.         {
  2735.             printf("Exp:%12f ", genBasis[i].exponents[j]);
  2736.             if (genBasis[i].s != NULL)
  2737.                 printf("S:%12f ", genBasis[i].s[j]);
  2738.             if (genBasis[i].p != NULL)
  2739.                 printf("P:%12f ", genBasis[i].p[j]);
  2740.             if (genBasis[i].d != NULL)
  2741.                 printf("D:%12f ", genBasis[i].d[j]);
  2742.             printf("\n");
  2743.         }
  2744.     }
  2745.     if ( normalize == true)
  2746.     {
  2747.         sum = 0.0;
  2748.         for ( i = 0; i < nBasis; i++ )
  2749.             sum += d[i] * d[i];
  2750.         N = 1.0 / sqrt(sum);        
  2751.         for ( i = 0; i < nBasis; i++ )
  2752.             d[i] = N * d[i];
  2753.         printf("Normalization coefficient: %f\n", N);
  2754.     } 
  2755.     printf("Molecular orbital coefficients for orbital #%d:\n", orbitalNumber);
  2756.     for ( i = 0; i < nBasis; i++)
  2757.         printf("%5d %10.5f\n", i+1, d[i]);
  2758.  
  2759.     /* ----------------------
  2760.         Now generate parameters 
  2761.         ---------------------- */
  2762.     
  2763.     if ( makeGSNParameters( ¶meters, atoms, nAtoms, nSym, nBasis, nPrim, \
  2764.     genBasis, nGen, shellsPerAtom, d) == -1 )
  2765.     {
  2766.         FileFormatError("ReadGSNFile", "Unable to generate parameters");
  2767.         return(NULLOBJ);
  2768.     }
  2769.     printf("GSN2 Parameters:\n");
  2770.     for ( i = 0; i < nPrim; i ++ )
  2771.         printf("%8.4f %8.4f %8.4f %5d %5d %5d %8.4f %8.4f %8.4f\n", parameters[i].d, parameters[i].c, \
  2772.         parameters[i].zeta, parameters[i].nx, parameters[i].ny, parameters[i].nz, \
  2773.         parameters[i].Rx,parameters[i].Ry, parameters[i].Rz );
  2774.     /* -----------------------------------
  2775.         Calculate the dimensions of the grid.
  2776.         There is always one more grid point than number
  2777.         of grid spaces in each direction.
  2778.        ------------------------------------ */
  2779.     dims[0] = ngridx + 1;
  2780.     dims[1] = ngridy + 1;
  2781.     dims[2] = ngridz + 1;
  2782.  
  2783.     /* -----------------------------------
  2784.         Create the scalar dataset
  2785.         ----------------------------------- */
  2786.     retVal = NewStructuredDataset(dataName, 3, dims, 0 );
  2787.  
  2788.     /* ------------------------------------
  2789.         Calculate grid bounds 
  2790.         ------------------------------------ */
  2791.     bounds[0] =  xmin;
  2792.     bounds[1] = xmax;
  2793.     bounds[2] =  ymin;
  2794.     bounds[3] = ymax;
  2795.     bounds[4] =  zmin;
  2796.     bounds[5] = zmax;
  2797.  
  2798.     /* ------------------------------------
  2799.         Create a regular data form
  2800.         ------------------------------------  */
  2801.     dataForm = NewRegularDataForm(formName, 3, dims, bounds);
  2802.     SetDatasetForm(retVal, dataForm);
  2803.     /* ------------------------------------
  2804.         Fill up the dataset (with x varying fastest)
  2805.         generating data values as you go
  2806.         Keep track of minimum and maximum values as
  2807.         data values are read in.
  2808.         ------------------------------------ */
  2809.     SetCurField(FIELD1, retVal);  /* Link the field to the data form */
  2810.     min = PLUSINF;                /* Initialize min and max to */
  2811.     max = MINUSINF;           /* extreme values */    
  2812.  
  2813.     x = xmin;
  2814.     y = ymin;
  2815.     z = zmin;
  2816.     xInc = (xmax - xmin ) / ( (real) (ngridx));
  2817.     yInc = (ymax - ymin) / ( (real) (ngridy));
  2818.     zInc = (zmax - zmin) / ((real) (ngridz));
  2819.  
  2820.     for (indices[2] = 0; indices[2] < dims[2]; (indices[2])++)
  2821.     {
  2822.         y = ymin;
  2823.         for (indices[1] = 0; indices[1] < dims[1]; (indices[1])++)
  2824.         {
  2825.             x = xmin;
  2826.             for (indices[0] = 0; indices[0] < dims[0]; (indices[0])++)
  2827.             {
  2828.                 curData = 0.0;
  2829.                 for (i=0; i < nPrim; i++ )
  2830.                     curData +=  parameters[i].d * parameters[i].c * \
  2831.                     gaussianFunction(x, y, z, parameters[i].Rx, parameters[i].Ry,\
  2832.                      parameters[i].Rz, parameters[i].zeta, parameters[i].nx, \
  2833.                      parameters[i].ny, parameters[i].nz, 1.0);
  2834.  
  2835.                 /*printf("%f ", curData);*/
  2836.  
  2837.                 PutFieldComponent(FIELD1, 0, indices,curData);
  2838.                 if (curData < min)
  2839.                     min = curData;
  2840.                 if (curData > max)
  2841.                     max = curData;
  2842.                 x += xInc;
  2843.             }
  2844.         /*    printf("\n");*/
  2845.             y += yInc;
  2846.         }
  2847.         z += zInc;
  2848.     }
  2849.  
  2850.     /* -------------------------
  2851.         Record the minimum and maximum data values
  2852.         for colour palette and isosurface scales 
  2853.         ------------------------- */
  2854.  
  2855.     minMax[0] = min;
  2856.     minMax[1] = max; 
  2857.     minMaxArray = NewRealArray(1, 2L);
  2858.     CArray2Array(minMaxArray, minMax);
  2859.     SetVar(retVal, MINMAX, minMaxArray);
  2860.  
  2861.     RegisterDataset(retVal);
  2862.     fclose(infile);
  2863.  
  2864.     /* -----------------------------
  2865.         Generate XYZ dataset
  2866.             ----------------------------- */
  2867.     makeXYZDataset( tmpName, atoms, nAtoms);
  2868.  
  2869.     return(NULLOBJ);
  2870.  
  2871.  
  2872. }
  2873. /* --------------------------------
  2874.     GSN File Reader Controls 
  2875.     Modified from AddNXRControls()
  2876.    -------------------------------- */
  2877. #define GSN_ITEMWIDTH 225
  2878. #define GSN_MIDDLE_SEP 45
  2879. #define GSN_SMALLER_ITEM 140
  2880. #define GSN_SMALLEST_ITEM 60
  2881. #define GSN_RADIOWIDTH 400
  2882. #define GSN_BOXLEFT 325
  2883. #define GSN_BOXWIDTH 60
  2884.  
  2885. #ifdef PROTO
  2886. static ObjPtr    AddGSNControls(ObjPtr FileReader, ObjPtr PanelContents)
  2887. #else
  2888. static ObjPtr    AddGSNControls(FileReader, PanelContents)
  2889.   ObjPtr  FileReader, PanelContents;
  2890. #endif
  2891. {
  2892.   ObjPtr titleBox, radio, checkBox, button;
  2893.   ObjPtr  checkbox, textbox, var;
  2894.   int     left, right, bottom, top;
  2895.  
  2896.   bottom = left = MAJORBORDER;
  2897.  
  2898. /* Grid bounds */
  2899. /* z */
  2900.   bottom += MINORBORDER;
  2901.   right = left + GSN_SMALLEST_ITEM;
  2902.   top = bottom + EDITBOXHEIGHT;
  2903.  
  2904.   textbox = NewTextBox(left, right, bottom, top,
  2905.                0, "Grid z min text", "z min:");
  2906.   SetVar(textbox, PARENT, PanelContents);
  2907.   PrefixList(PanelContents, textbox);
  2908.  
  2909.     left = right;
  2910.   right = left + GSN_BOXWIDTH;
  2911.   var = GetRealVar("AddGSNControls", FileReader, GSNZMIN);
  2912.   if (!var)
  2913.   {
  2914.     SetVar(FileReader, GSNZMIN, NewReal(GZMIN));
  2915.  
  2916.   }
  2917.  
  2918.  
  2919.   textbox = NewTextBox(left, right, bottom, top,
  2920.                EDITABLE + WITH_PIT + ONE_LINE,
  2921.                "zmin", "0");
  2922.   SetVar(textbox, PARENT, PanelContents);
  2923.   PrefixList(PanelContents, textbox);
  2924.   SetTextAlign(textbox, RIGHTALIGN);
  2925.   AssocTextRealControlWithVar(textbox, FileReader, GSNZMIN, 
  2926.                   minusInf, plusInf,TR_NE_TOP);
  2927.   SetVar(textbox, HELPSTRING,
  2928.      NewString("This value is the minimum z value of grid bounds."));
  2929.  
  2930.   left = right + GSN_MIDDLE_SEP;
  2931.   right = left + GSN_SMALLEST_ITEM;
  2932.   textbox = NewTextBox(left, right, bottom, top,
  2933.                0, "Grid z max text", "z max:");
  2934.   SetVar(textbox, PARENT, PanelContents);
  2935.   PrefixList(PanelContents, textbox);
  2936.  
  2937.     left = right;
  2938.   right = left + GSN_BOXWIDTH;
  2939.   var = GetRealVar("AddGSNControls", FileReader, GSNZMAX);
  2940.   if (!var)
  2941.     SetVar(FileReader, GSNZMAX, NewReal(GZMAX));
  2942.   
  2943.   textbox = NewTextBox(left, right, bottom, top,
  2944.                EDITABLE + WITH_PIT + ONE_LINE,
  2945.                "zmax", "0");
  2946.   SetVar(textbox, PARENT, PanelContents);
  2947.   PrefixList(PanelContents, textbox);
  2948.   SetTextAlign(textbox, RIGHTALIGN);
  2949.  
  2950.   AssocTextRealControlWithVar(textbox, FileReader, GSNZMAX, 
  2951.                  minusInf , plusInf,TR_NE_TOP);
  2952.   SetVar(textbox, HELPSTRING,
  2953.      NewString("This value is the maximum z value of grid bounds."));
  2954.  
  2955.   left = right + GSN_MIDDLE_SEP;
  2956.   right = left + GSN_SMALLER_ITEM;
  2957.   textbox = NewTextBox(left, right, bottom, top,
  2958.                0, "Grid z spacing", "Grid spacing(z):");
  2959.   SetVar(textbox, PARENT, PanelContents);
  2960.   PrefixList(PanelContents, textbox);
  2961.  
  2962.     left = right;
  2963.   right = left + GSN_BOXWIDTH;
  2964.   var = GetIntVar("AddGSNControls", FileReader, GSNGRIDZ);
  2965.   if (!var)
  2966.     SetVar(FileReader, GSNGRIDZ, NewInt(GGRIDZ));
  2967.   textbox = NewTextBox(left, right, bottom, top,
  2968.                EDITABLE + WITH_PIT + ONE_LINE,
  2969.                "gridz", "0");
  2970.   SetVar(textbox, PARENT, PanelContents);
  2971.   PrefixList(PanelContents, textbox);
  2972.   SetTextAlign(textbox, RIGHTALIGN);
  2973.   AssocTextIntControlWithVar(textbox, FileReader, GSNGRIDZ, 
  2974.                   0.0, plusInf, TR_INT_ONLY +TR_NE_TOP);
  2975.   SetVar(textbox, HELPSTRING,
  2976.      NewString("This value is the grid spacing in the z direction."));
  2977.   
  2978.   bottom = top + MINORBORDER;
  2979.   left = MAJORBORDER;
  2980.   top = bottom + EDITBOXHEIGHT;
  2981.   right = left + GSN_SMALLEST_ITEM;
  2982.  
  2983.   textbox = NewTextBox(left, right, bottom, top,
  2984.                0, "Grid y min text", "y min:");
  2985.   SetVar(textbox, PARENT, PanelContents);
  2986.   PrefixList(PanelContents, textbox);
  2987.  
  2988.     left = right;
  2989.   right = left + GSN_BOXWIDTH;
  2990.   var = GetRealVar("AddGSNControls", FileReader, GSNYMIN);
  2991.   if (!var)
  2992.     SetVar(FileReader, GSNYMIN, NewReal(GYMIN));
  2993.   textbox = NewTextBox(left, right, bottom, top,
  2994.                EDITABLE + WITH_PIT + ONE_LINE,
  2995.                "ymin", "0");
  2996.   SetVar(textbox, PARENT, PanelContents);
  2997.   PrefixList(PanelContents, textbox);
  2998.   SetTextAlign(textbox, RIGHTALIGN);
  2999.   AssocTextRealControlWithVar(textbox, FileReader, GSNYMIN, 
  3000.                   minusInf, plusInf,TR_NE_TOP);
  3001.   SetVar(textbox, HELPSTRING,
  3002.      NewString("This value is the minimum y value of grid bounds."));
  3003.  
  3004.   left = right + GSN_MIDDLE_SEP;
  3005.   right = left + GSN_SMALLEST_ITEM;
  3006.   textbox = NewTextBox(left, right, bottom, top,
  3007.                0, "Grid y max text", "y max:");
  3008.   SetVar(textbox, PARENT, PanelContents);
  3009.   PrefixList(PanelContents, textbox);
  3010.  
  3011.     left = right;
  3012.   right = left + GSN_BOXWIDTH;
  3013.   var = GetRealVar("AddGSNControls", FileReader, GSNYMAX);
  3014.   if (!var)
  3015.     SetVar(FileReader, GSNYMAX, NewReal(GYMAX));
  3016.   textbox = NewTextBox(left, right, bottom, top,
  3017.                EDITABLE + WITH_PIT + ONE_LINE,
  3018.                "ymax", "0");
  3019.   SetVar(textbox, PARENT, PanelContents);
  3020.   PrefixList(PanelContents, textbox);
  3021.   SetTextAlign(textbox, RIGHTALIGN);
  3022.   AssocTextRealControlWithVar(textbox, FileReader, GSNYMAX, 
  3023.                   minusInf, plusInf,TR_NE_TOP);
  3024.   SetVar(textbox, HELPSTRING,
  3025.      NewString("This value is the maximum y value of grid bounds."));
  3026.  
  3027.   left = right + GSN_MIDDLE_SEP;
  3028.   right = left + GSN_SMALLER_ITEM;
  3029.   textbox = NewTextBox(left, right, bottom, top,
  3030.                0, "Grid y spacing", "Grid spacing(y):");
  3031.   SetVar(textbox, PARENT, PanelContents);
  3032.   PrefixList(PanelContents, textbox);
  3033.  
  3034.     left = right;
  3035.   right = left + GSN_BOXWIDTH;
  3036.   var = GetIntVar("AddGSNControls", FileReader, GSNGRIDY);
  3037.   if (!var)
  3038.     SetVar(FileReader, GSNGRIDY, NewInt(GGRIDY));
  3039.   textbox = NewTextBox(left, right, bottom, top,
  3040.                EDITABLE + WITH_PIT + ONE_LINE,
  3041.                "gridy", "0");
  3042.   SetVar(textbox, PARENT, PanelContents);
  3043.   PrefixList(PanelContents, textbox);
  3044.   SetTextAlign(textbox, RIGHTALIGN);
  3045.   AssocTextIntControlWithVar(textbox, FileReader, GSNGRIDY, 
  3046.                   0.0, plusInf, TR_INT_ONLY +TR_NE_TOP);
  3047.   SetVar(textbox, HELPSTRING,
  3048.      NewString("This value is the grid spacing in the y direction."));
  3049.  
  3050. /* x */
  3051.   left = MAJORBORDER;    
  3052.   bottom = top + MINORBORDER;
  3053.   right = left + GSN_SMALLEST_ITEM;
  3054.   top = bottom + EDITBOXHEIGHT;
  3055.  
  3056.   textbox = NewTextBox(left, right, bottom, top,
  3057.                0, "Grid x min text", "x min:");
  3058.   SetVar(textbox, PARENT, PanelContents);
  3059.   PrefixList(PanelContents, textbox);
  3060.  
  3061.     left = right;
  3062.   right = left + GSN_BOXWIDTH;
  3063.   var = GetRealVar("AddGSNControls", FileReader, GSNXMIN);
  3064.   if (!var)
  3065.     SetVar(FileReader, GSNXMIN, NewReal(GXMIN));
  3066.   textbox = NewTextBox(left, right, bottom, top,
  3067.                EDITABLE + WITH_PIT + ONE_LINE,
  3068.                "xmin", "0");
  3069.   SetVar(textbox, PARENT, PanelContents);
  3070.   PrefixList(PanelContents, textbox);
  3071.   SetTextAlign(textbox, RIGHTALIGN);
  3072.   AssocTextRealControlWithVar(textbox, FileReader, GSNXMIN, 
  3073.                   minusInf, plusInf,TR_NE_TOP);
  3074.   SetVar(textbox, HELPSTRING,
  3075.      NewString("This value is the minimum x value of grid bounds."));
  3076.  
  3077.   left = right + GSN_MIDDLE_SEP;
  3078.   right = left + GSN_SMALLEST_ITEM;
  3079.   textbox = NewTextBox(left, right, bottom, top,
  3080.                0, "Grid x max text", "x max:");
  3081.   SetVar(textbox, PARENT, PanelContents);
  3082.   PrefixList(PanelContents, textbox);
  3083.  
  3084.     left = right;
  3085.   right = left + GSN_BOXWIDTH;
  3086.   var = GetRealVar("AddGSNControls", FileReader, GSNXMAX);
  3087.   if (!var)
  3088.     SetVar(FileReader, GSNXMAX, NewReal(GXMAX));
  3089.   textbox = NewTextBox(left, right, bottom, top,
  3090.                EDITABLE + WITH_PIT + ONE_LINE,
  3091.                "xmax", "0");
  3092.   SetVar(textbox, PARENT, PanelContents);
  3093.   PrefixList(PanelContents, textbox);
  3094.   SetTextAlign(textbox, RIGHTALIGN);
  3095.   AssocTextRealControlWithVar(textbox, FileReader, GSNXMAX, 
  3096.                   minusInf, plusInf,TR_NE_TOP);
  3097.   SetVar(textbox, HELPSTRING,
  3098.      NewString("This value is the maximum x value of grid bounds."));
  3099.  
  3100.  
  3101.   left = right + GSN_MIDDLE_SEP;
  3102.   right = left + GSN_SMALLER_ITEM;
  3103.   textbox = NewTextBox(left, right, bottom, top,
  3104.                0, "Grid x spacing", "Grid spacing(x):");
  3105.   SetVar(textbox, PARENT, PanelContents);
  3106.   PrefixList(PanelContents, textbox);
  3107.  
  3108.     left = right;
  3109.   right = left + GSN_BOXWIDTH;
  3110.   var = GetIntVar("AddGSNControls", FileReader, GSNGRIDX);
  3111.   if (!var)
  3112.     SetVar(FileReader, GSNGRIDX, NewInt(GGRIDX));
  3113.   textbox = NewTextBox(left, right, bottom, top,
  3114.                EDITABLE + WITH_PIT + ONE_LINE,
  3115.                "gridx", "0");
  3116.   SetVar(textbox, PARENT, PanelContents);
  3117.   PrefixList(PanelContents, textbox);
  3118.   SetTextAlign(textbox, RIGHTALIGN);
  3119.   AssocTextIntControlWithVar(textbox, FileReader, GSNGRIDX, 
  3120.                   0.0, plusInf, TR_INT_ONLY +TR_NE_TOP);
  3121.   SetVar(textbox, HELPSTRING,
  3122.      NewString("This value is the grid spacing in the x direction."));
  3123.  
  3124.  
  3125.   left = MAJORBORDER;
  3126.   bottom = top + EDITBOXHEIGHT;
  3127.  
  3128. /* Normalization of coefficients checkbox */
  3129.     top = bottom + CHECKBOXHEIGHT;
  3130.     right = left + FRTIMEWIDTH;
  3131.     checkBox = NewCheckBox(left, right, bottom, top,
  3132.         "Normalize coefficients",
  3133.         GetPredicate(FileReader, USER8));
  3134.     SetVar(checkBox, PARENT, PanelContents);
  3135.     PrefixList(PanelContents, checkBox);
  3136.     bottom = top + MINORBORDER;
  3137.     if (!GetVar(FileReader, USER8))
  3138.     {
  3139.         SetVar(FileReader, USER8, ObjFalse);
  3140.     }
  3141.     AssocDirectControlWithVar(checkBox, FileReader, USER8);
  3142.     SetVar(checkBox, HELPSTRING, 
  3143.         NewString("If this box is checked, the coefficients for the \
  3144. molecular orbital selected will be normalized such that the sum of \
  3145. the squares of the coefficients equals 1."));
  3146.  
  3147.  
  3148.  /* Molecular orbital textbox*/
  3149.   right = left + GSN_SMALLER_ITEM;
  3150.   top = bottom + EDITBOXHEIGHT;
  3151.  
  3152.   textbox = NewTextBox(left, right, bottom, top,
  3153.                0, "Molecular Orbital Text", "Molecular Orbital #:");
  3154.   SetVar(textbox, PARENT, PanelContents);
  3155.   PrefixList(PanelContents, textbox);
  3156.  
  3157.  
  3158.  
  3159.     left = right;
  3160.   right = left + GSN_BOXWIDTH;
  3161.  
  3162.  
  3163.   var = GetIntVar("AddGSNControls", FileReader, USER7);
  3164.   if (!var)
  3165.     SetVar(FileReader, USER7, NewReal(MOLECULARORBITAL));
  3166.   textbox = NewTextBox(left, right, bottom, top,
  3167.                EDITABLE + WITH_PIT + ONE_LINE,
  3168.                "Molecular Orbital #", "0");
  3169.   SetVar(textbox, PARENT, PanelContents);
  3170.   PrefixList(PanelContents, textbox);
  3171.   SetTextAlign(textbox, RIGHTALIGN);
  3172.   AssocTextIntControlWithVar(textbox, FileReader, USER7, 
  3173.                   1.0, plusInf, TR_INT_ONLY + TR_NE_TOP);
  3174.   SetVar(textbox, HELPSTRING,
  3175.      NewString("This value is the molecular orbital # selected."));
  3176.  
  3177.  
  3178.   return NULLOBJ;
  3179. }
  3180.  
  3181. /* ------------------------------------
  3182.     Initialize SIMS File Readers
  3183.     ------------------------------------ */
  3184. #ifdef PROTO
  3185. void    InitSIMSFiles(void)
  3186. #else
  3187. void    InitSIMSFiles()
  3188. #endif
  3189. {
  3190.  
  3191.   ObjPtr  fileReader;
  3192.  
  3193.   fileReader = NewFileReader("XYZ");
  3194.   SetVar(fileReader, EXTENSION, NewString("xyz"));
  3195.   SetVar(fileReader, HELPSTRING, \
  3196.     NewString("This file reader takes molecular data in xyz format as specified by\
  3197.  the Minnesota Supercomputer Center's XMol.  Format specification follows: \
  3198. There may be multiple frames for animation.  For each frame, the first line \
  3199. must have an integer specifying the number of atoms, the second line is reserved \
  3200. for a title (possibly just a blank line), and one line for each atom in the form: \
  3201. name x y z (charge) (vectorX) (vectorY) (vectorZ).  A number of frames may \
  3202. follow. Charge and/or field vector \
  3203. may be specified but are ignored by the reader.  It is assumed that the number of \
  3204. atoms will be the same from frame to frame.  The file reader produces three \
  3205. datasets, with the extensions: .atoms, .radii, .form.  The radii and form datasets \
  3206. may be ignored by the user.  The atoms dataset may be visualized in ball and or \
  3207. stick form.  Developed \
  3208. at the Steacie Institute for Molecular Sciences at the National Research Council \
  3209. of Canada."));
  3210.   SetMethod(fileReader, READALL, ReadXYZFile);
  3211.   SetMethod(fileReader, ADDCONTROLS, AddXYZControls);
  3212.  
  3213. /* -------------------------------------
  3214.     set USER1 to be the minimum distance required for a bond
  3215.    ------------------------------------- */
  3216.   SetVar(fileReader, USER1, NewReal(MAXSTICKLENGTH)); 
  3217.  
  3218. /* --------------------------------------
  3219.     set USER3 to be radio button group for bond calculation
  3220.     -------------------------------------- */
  3221.   SetVar(fileReader, USER3, NewInt(CALC_COVALENTRADII));
  3222. /* ---------------------------------------
  3223.     set USER4 to be the scaling factor used in covalent radii calculation
  3224.     --------------------------------------- */
  3225.   SetVar(fileReader, USER4, NewReal(FUDGEFACTOR));
  3226. /* ---------------------------------------
  3227.     set USER5 to be pre-set scaling factor of sphere radius
  3228.     e.g. want to set this low for reading in a lot of organic molecules
  3229.     --------------------------------------- */
  3230.   SetVar(fileReader, USER5, NewReal(1.0));
  3231.  
  3232. /*----------------------------------------
  3233.     set USER6 to be checkbox for conversion from atomic
  3234.     length units to angstroms
  3235.    ---------------------------------------- */
  3236.   SetVar(fileReader, USER6, ObjFalse);
  3237.  
  3238.  
  3239.   /* save setting stuff */
  3240.   AddSnapVar(fileReader, USER1);
  3241.   AddSnapVar(fileReader, USER2);
  3242.   AddSnapVar(fileReader, USER3);
  3243.   AddSnapVar(fileReader, USER4);
  3244.   AddSnapVar(fileReader, USER5);
  3245.   AddSnapVar(fileReader, USER6);
  3246.  
  3247.   ApplySavedSettings(fileReader);
  3248.  
  3249. /* -----------------------------------------
  3250.     Read in user-defined atoms if any exist
  3251.     ----------------------------------------- */
  3252.    numExtraXYZSymbols = ReadExtraXYZSymbolTable();
  3253.    if (numExtraXYZSymbols== -1)
  3254.     FileFormatError("ReadExtraXYZSymbolTable", "Poor format for user defined XYZ \
  3255. symbol table");
  3256.  
  3257.  
  3258.   fileReader = NewFileReader("GRD");
  3259.   SetVar(fileReader, EXTENSION, NewString("grd"));
  3260.   SetVar(fileReader, HELPSTRING, \
  3261.     NewString("This file reader takes 3D grid data in .grd format as generated by \
  3262. BioSym's DMol PLOT command (see p.B31-B35 of DMol user manual).  \
  3263. Developed at the Steacie Institute for Molecular Sciences at the \
  3264. National Research Council of Canada."));
  3265.   SetMethod(fileReader, READALL, ReadGRDFile);
  3266.  
  3267.   fileReader = NewFileReader("GSN1");
  3268.   SetVar(fileReader, EXTENSION, NewString("gsn1"));
  3269.   SetVar(fileReader, HELPSTRING, \
  3270.     NewString("The purpose of the GSN1 file reader is to take \
  3271. as input the coefficients and exponents of gaussian primitives, \
  3272. obtained from ab initio calculations such as that produced \
  3273. by G92, and approximate the molecular orbitals by a linear \
  3274. combination of the primitives. \n\n \
  3275. The value computed at each grid point is obtained from \
  3276. the function:\n\n  f(x,y,z) = SUM(i=1..n){ c_i * g_i (zeta,x,y,z,nx,ny,nz, \
  3277. Rx, Ry, Rz ) }\n \
  3278.  The GSN1 file must be of the form:\n \
  3279. **************************************************\n \
  3280. TITLE\n \
  3281. xmin xmax ymin ymax zmin zmax\n \
  3282. nGridx nGridy nGridz\n \
  3283. n\n \
  3284. c1 zeta1 nx1 ny1 nz1 Rx1 Ry1 Rz1\n \
  3285. c2 zeta2 nx2 ny2 nz2 Rx2 Ry2 Rz2\n \
  3286. .\n.\n.\n \
  3287. cn zetan nxn nyn nzn Rxn Ryn Rzn\n \
  3288. *************************************************\n \
  3289. Developed at the Steacie Institute for Molecular Sciences at the \
  3290. National Research Council of Canada."));
  3291.  
  3292.  
  3293.   SetMethod(fileReader, READALL, ReadGSN1File);
  3294.  
  3295.   fileReader = NewFileReader("GSN2");
  3296.   SetVar(fileReader, EXTENSION, NewString("gsn2"));
  3297.   SetVar(fileReader, HELPSTRING, \
  3298.     NewString("The purpose of the GSN2 file reader is to take \
  3299. as input the coefficients and exponents of gaussian primitives, \
  3300. obtained from ab initio calculations such as that produced \
  3301. by G92, and approximate the molecular orbitals by a linear \
  3302. combination of the primitives. \n\n \
  3303. The value computed at each grid point is obtained from \
  3304. the function:\n\n  f(x,y,z) = SUM(i=1..n){ d_i * c_i * g_i (zeta,x,y,z,nx,ny,nz, \
  3305. Rx, Ry, Rz ) }\n \
  3306.  The GSN2 file must be of the form:\n \
  3307. **************************************************\n \
  3308. TITLE\n \
  3309. xmin xmax ymin ymax zmin zmax\n \
  3310. nGridx nGridy nGridz\n \
  3311. n\n \
  3312. d1 c1 zeta1 nx1 ny1 nz1 Rx1 Ry1 Rz1\n \
  3313. d2 c2 zeta2 nx2 ny2 nz2 Rx2 Ry2 Rz2\n \
  3314. .\n.\n.\n \
  3315. dn cn zetan nxn nyn nzn Rxn Ryn Rzn\n \
  3316. *************************************************\n \
  3317. Developed at the Steacie Institute for Molecular Sciences at the \
  3318. National Research Council of Canada."));
  3319.  
  3320.  SetMethod(fileReader, READALL, ReadGSN2File);
  3321.  
  3322.   fileReader = NewFileReader("GSN");
  3323.   SetVar(fileReader, EXTENSION, NewString("gsn"));
  3324.   SetVar(fileReader, HELPSTRING, \
  3325.     NewString("The GSN file reader will extract the coefficients \
  3326. and exponents of gaussian primitives from a G92 output file, and \
  3327. approximate the molecular orbitals by a linear combination of the \
  3328. primitives. In addition, it will take the atomic coordinate information \
  3329. in the G92 output file to generate a .atoms dataset for visualizing \
  3330. the molecule in ball-and-or-stick form.  The combined visualization \
  3331. of the molecular orbitals and the ball/stick molecule can be \
  3332. quite informative. \n\nGaussian92 File Specifications\n \
  3333. ===============================\n \
  3334. The GSN file reader can read standard N-M1G basis sets, with\
  3335. or without additional polarization or diffuse functions.\
  3336. The reader will recognize S,P,D shell types, but not F orbitals.\
  3337. \n \
  3338. To create an output file that may be read by the GSN reader,\
  3339. the command line of the Gaussian92 input must contain\
  3340. \na) the keyword to print the GAUSSIAN function table in \
  3341. GenBas format:\niop(3/24=10)\n \
  3342. b) the keyword to print eigenvalues and eigenvectors:\n \
  3343. Pop=reg {or Pop=full}\n \
  3344. Developed at the Steacie Institute for Molecular Sciences at the \
  3345. National Research Council of Canada."));
  3346.  
  3347.   SetMethod(fileReader, READALL, ReadGSNFile);
  3348.   SetMethod(fileReader, ADDCONTROLS, AddGSNControls);
  3349.  
  3350. /* -------------------------------------
  3351.     set USER7 to be the molecular orbital selected
  3352.     set USER8 to be checkbox for normalization of
  3353.     molecular orbital coefficients (default: off)
  3354.    ------------------------------------- */
  3355.   SetVar(fileReader, USER7, NewInt(MOLECULARORBITAL)); 
  3356.   SetVar(fileReader, USER8, ObjFalse);
  3357.  
  3358. /* --------------------------------------
  3359.     grid bounds/spacing variables
  3360.    -------------------------------------- */
  3361.   SetVar(fileReader, GSNXMIN, NewReal(GXMIN));
  3362.   SetVar(fileReader, GSNXMAX, NewReal(GXMAX));
  3363.   SetVar(fileReader, GSNYMIN, NewReal(GYMIN));
  3364.   SetVar(fileReader, GSNYMAX, NewReal(GYMAX));
  3365.   SetVar(fileReader, GSNZMIN, NewReal(GZMIN));
  3366.   SetVar(fileReader, GSNZMAX, NewReal(GZMAX));
  3367.   
  3368.   SetVar(fileReader, GSNGRIDX, NewInt(GGRIDX));
  3369.   SetVar(fileReader, GSNGRIDY, NewInt(GGRIDY));
  3370.   SetVar(fileReader, GSNGRIDZ, NewInt(GGRIDZ));
  3371.   
  3372.  /* save setting stuff */
  3373.   AddSnapVar(fileReader, USER7);
  3374.   AddSnapVar(fileReader, USER8);
  3375.   AddSnapVar(fileReader, GSNXMIN);
  3376.   AddSnapVar(fileReader, GSNXMAX);
  3377.   AddSnapVar(fileReader, GSNYMIN);
  3378.   AddSnapVar(fileReader, GSNYMAX);
  3379.   AddSnapVar(fileReader, GSNZMIN);
  3380.   AddSnapVar(fileReader, GSNZMAX);
  3381.   AddSnapVar(fileReader, GSNGRIDX);
  3382.   AddSnapVar(fileReader, GSNGRIDY);
  3383.   AddSnapVar(fileReader, GSNGRIDZ);
  3384.  
  3385.   ApplySavedSettings(fileReader);
  3386.  
  3387.  
  3388. }
  3389.  
  3390. /* ------------------------------------
  3391.     Kills any SIMS File Handling Routines
  3392.     ------------------------------------ */
  3393. void KillSIMSFiles()
  3394. {
  3395.  
  3396. }
  3397.